State-Dependent Molecular Dynamics

This paper proposes a new mixed quantum mechanics (QM)—molecular mechanics (MM) approach, where MM is replaced by quantum Hamilton mechanics (QHM), which inherits the modeling capability of MM, while preserving the state-dependent nature of QM. QHM, a single mechanics playing the roles of QM and MM simultaneously, will be employed here to derive the three-dimensional quantum dynamics of diatomic molecules. The resulting state-dependent molecular dynamics including vibration, rotation and spin are shown to completely agree with the QM description and well match the experimental vibration-rotation spectrum. QHM can be incorporated into the framework of a mixed quantum-classical Bohmian method to enable a trajectory interpretation of orbital-spin interaction and spin entanglement in molecular dynamics.


Introduction
One of the greatest challenges in molecular dynamics (MD) is to model processes involving many degrees of freedom, some of which have to be treated quantum mechanically. The combined quantum mechanics (QM) and molecular mechanics (MM) approach to MD [1,2] provides tremendous computational advantages over full quantum mechanical models by treating a limited region of a molecular system quantum mechanically, while treating the rest of the system by using conventional OPEN ACCESS MM methods. The QM part corresponds to what is to be studied in detail, such as the regions involving charge transfer, electron excitations or chemical reactions. The atoms in the QM part are explicitly expressed as electrons and nuclei such that their motions are described quantum mechanically. Quantum effects in MD simulation have been taken into account by various techniques, such as pseudo-spectral methods [3], path integral methods [4], Car-Parrinello-type simulations [5], molecular wave packet dynamics [6] and density matrix evolution method [7]. The performance of a hybrid QM/MM method has been extensively investigated by Billing and coworkers [8]. Brickmann and Schmitt [9] employed a sequence of approximations to transfer the mixed QM/MM dynamics into a Hamilton-Jacobi-type scheme, which is able to formulate the equations of motion by using a single Hamiltonian.
The QM/MM approach to MD simulation of molecular systems mainly relies on the Born-Oppenheimer approximation that a molecular motion can be separated into two independent motions: one is the classical atomic motion on a single adiabatic potential energy surface and the other is the quantum electronic motion in the presence of a time-dependent potential generated by the moving classical atoms. However, there are a huge number of problems for which the interaction between classical and quantum motion is so significant that the Born-Oppenheimer approximation may become invalid. For such problems, the QM/MM approach to MD simulation inevitably faces the crucial issue of self-consistency. The quantum degrees of freedom must evolve correctly under the influence of the surrounding classical motion, and meanwhile, the classical degrees of freedom must respond correctly to the quantum transitions. Therefore, a self-consistent QM/MM approach has to consider not only the effect of the classical degrees of freedom on the quantum ones but also the backreaction of the quantum effect on the classical degrees of freedom.
Two mixed quantum-classical methods that have been developed in the literature to treat the interactions between classical and quantum systems in a self-consistent way are the mean-field method [10] and the surface-hopping method [11]. The mean-field method calculates the force for the classical motion by averaging over the quantum wave function. This method is invariant to the choice of quantum representations and applicable to both bound and continuum states, but it suffers from the neglect of correlations between classical and quantum degree of freedom. On the other hand, the surface-hopping method was developed to manifest quantum-classical correlation, but it is not invariant to the choice of quantum representations and is intrinsically limited to discrete quantum states. Although both methods have their respective limitations, they have been proven to be successful in many applications.
In the latest decade, a novel solution to the quantum backreaction problem in a mixed QM/MM simulation has been proposed using the Bohmian formulation of quantum mechanics [12,13]. The mixed quantum-classical Bohmian (MQCB) approach combines the merits of the above two methods and gives a consistent treatment of mixing quantum and classical degrees of freedom without reference to any basis set. The MQCB method has been applied to the process of vibrational decoherence of I2 in a dense helium environment [14], to the case of rotational diffractive surface scattering of a diatomic molecule [15] and to a model of O2 interacting with a Pt surface [16], all with good agreement with the full quantum-mechanical treatments.
The current formulation of Bohmian mechanic is mainly based on Cartesian coordinates, which is inconvenient to the description of molecular angular motions. This paper aims to incorporate quantum Hamilton mechanics (QHM) [17][18][19] into the framework of the MQCB method to enhance its capability of handling orbital and spin angular motions. It is well known that the Bohmian velocity vanishes in all of the stationary states with zero angular quantum number, and therefore, the corresponding Bohmian particles will remain standing at their initial positions at any time [20]. QHM, on the one hand, conquers the problem of Bohmian stationarity by formulating Bohmian mechanics in complex space, and on the other, describes MD in terms of Hamilton equations of motion, which are coordinate-independent and suitable to any curvilinear coordinates.
With the proposed modifications by QHM, the computational procedures of MQCB developed in the literature can be used to simulate molecular dynamics, including vibration, rotation and spin motions. The correctness of the derived state-dependent molecular dynamics will be verified by comparing with the quantum mechanical description of a diatomic molecule for which the Schrödinger equation has an analytical solution. Of significance is that the spin dynamics, which cannot be described by spatial wave functions, emerges naturally from the established state-dependent molecular dynamics.
In Section 2, we introduce the working equations of the MQCB method to describe the motion of a diatomic molecule in Cartesian coordinates. QHM is then introduced in Section 3 to reformulate Bohmian mechanics under spherical coordinates. We show that QHM is a single mechanics simultaneously playing the roles of QM and MM by pointing out that QHM comprises two sets of Hamilton equations with the first set describing the vibration-rotation motion of the molecule and the second set describing the time evolution of the wave function. State-dependent vibrational dynamics fully consistent with QM is derived from QHM in Section 4. Section 5 demonstrates how to incorporate molecular spin dynamics into MD simulation under the same motion space governed by QHM. The resulting state-dependent molecular dynamics is found to agree with the prediction of QM and well match the experimental vibration-rotation spectrum.

Mixed Quantum-Classical Mechanics
In this section, we review the basics of the MQCB method by applying it to simulate the motion of a diatomic molecule, which is composed of Nuclei A and B, of mass MA and MB together with a number N of electrons (see Figure 1). The internuclear position vector will be denoted by X and the position vectors of the electrons with respect to O, the center of mass of A and B, by r 1 , r 2 ,⋯, r N . The position 2) can be recast into classical-like equations of motion: A workable scheme of MQCB has been implemented successfully by replacing the full-dimensional wave function ( , , ) by a wave function ( , , ) in the quantum subspace that depends parametrically on the classical position ( ). The approximate wave function ( , , ) obeys the Schrödinger equation: 2 2 2 ( , , ( )) ( , ( )) ( , , ( )) 2 where / denotes the material derivative / + • along the classical trajectory ( ) . The MQCB method and the mean-field method share the same quantum degree of freedom described by Equation (2.4), but differ from each other in the way that the classical degree of freedom is governed by the quantum wave function. Instead of Equation (2.5), the classical trajectory ( ) for the mean-field method is computed by the mean potential: Upon applying Equation (2.5) or Equation (2.7) to a diatomic molecule, we find that it is convenient to express the vector in the spherical coordinates ( , , ) in order to describe the vibration and rotation of the two atoms with respect to their center of mass. However, the expressions of quantum motion described by Equation (2.2) to Equation (2.7) are valid only for Cartesian coordinates. To give an explicit manifestation of the orbital and spin angular motion of the molecules, the reformulation of Bohmian mechanics under spherical coordinates is necessary. It will be shown in the next section that QHM provides MQCB with state-dependent molecular dynamics, which preserve the quantization property of orbital and spin angular momentum. By contrast, in the mean-field method, the classical trajectories ( ) over different quantum states are averaged out to give a state-independent result, as can be seen from Equation (2.7). To compare with the QHM formulation discussed later, here let us have a quick look on the description of the vibration-rotation motion of a diatomic molecule by the classical Equation (2.7). With the effective potential ( , ) = ( , , , ) given by Equation (2.7), the classical Hamiltonian governing the 3D relative motion of the two nuclei with respect to their center of mass can be expressed by: . (2.8) The molecular dynamics ( ( ), ( ), ( )) can be solved from the Hamilton equations: , It is clear that the vibrational and rotational motions of the molecule are determined uniquely by the effective potential ( , , , ) with initial conditions ( (0), (0), (0)) and ( (0), (0), (0)).
In the case of a central-force potential ( ), the total energy , the squared angular momentum and the z-component angular momentum are conserved quantities along any dynamic trajectory ( ( ), ( ), ( )): , (2.11) where the three conservation constants depend continuously on the initial conditions. In the next section, a quantum version of the classical Hamilton Equations (2.9) and (2.10) will be derived by QHM. The resulting quantum Hamilton equations can describe molecular dynamics in a general non-stationary quantum state Ψ( , , , ), and a quantum version of the conservation law (2.11) comes out naturally, if Ψ( , , , ) is an eigenfunction.

Quantum Hamilton Mechanics
Consider the diatomic molecule as shown in Figure 1b. Let be the internuclear position vector expressed in a curvilinear coordinate system, ( , ) be the instantaneous internuclear potential, and be the reduced mass. The Schrödinger equation describing the molecular motion can be recast into the quantum Hamilton-Jacobi equation [17,20]: where the action function is the logarithmic wave function defined as: . (3. 2) The canonical momentum in Equation (3.1) is related to the action function via the law of canonical transformation: .
The appearance of the imaginary number = √−1 indicates that the momentum has to be defined in a complex domain. In Cartesian coordinates, we have = M and the particle's velocity can be determined by the wave function Ψ as , (3.4) which is the governing equation in the complex-valued Bohmian mechanics [21]. It is the complex momentum from Equation (3.3) rather than the real momentum from Equation (2.2) that matches the momentum distribution provided by standard quantum mechanics [20]. The relationship between the real and complex momentum can be found as: It can be seen that the real Bohmian trajectory solved from = (1/ ) / only carry information about the dynamics of the momentum flow, while the complex trajectory solved from = (1/ ) / also includes information about the probability = |Ψ|. The dynamics in the complex phase space ( , ) thus explains in a natural way how to get the correct momentum distribution and explains why algorithms based on complex trajectories are stable and accurate [22][23][24][25][26][27][28]. The advantages of implementing numerical codes in a complex phase space are similar to those of considering complex fields instead of real ones in electromagnetism [20].
When compared with the classical Hamiltonian in Equation (2.8), the quantum Hamiltonian defined in Equation (3.1) has an additional term called complex quantum potential, which in the state Ψ can be expressed by: . (3.6) The state-dependent molecular dynamics to be developed here all originate from the state-dependent nature of . With the quantum Hamiltonian defined in Equation (3.1), the accompanying Hamilton equations appear as the usual form: The Hamilton equations are coordinate-independent and valid for arbitrary coordinate systems. Especially, in the Cartesian coordinates, we have ∂ / = / and the first set of the Hamilton equations (3.7) turns out to be the Bohmian velocity, like the one defined in Equation (2.2). We now specialize to the spherical coordinates = ( , , ) with denoting the internuclear distance, and ( , ) denoting the orientation of the molecule. By expanding the inner product • and the divergence ∇ • in spherical coordinates, quantum Hamiltonian in Equation (3.1) can be expressed by [17]: This is the quantum counterpart of the classical Hamiltonian defined in Equation (2.8). Substituting into Equation (3.7), we obtain the first set of the Hamilton equations as: where the canonical momentum ( , , ) has been given according to Equation , where the first three terms constitute the classical kinetic energy , and the last two terms form the total potential: The dynamic equations of motion for the molecule can be described either by the Hamilton equations (3.9) or by the Lagrange equations based on the quantum Lagrangian: , (3.13) from which the quantum Lagrange equations of motion can be obtained as It can be seen that apart from the internuclear force produced by , there are additional quantum forces acting on the molecule yielded by the quantum potential . In the next two sections, we proceed to show that the action of the quantum potential leads to the state-dependent molecular dynamics compatible with the description of the wave function Ψ.

State-Dependent Molecular Vibration
By replacing Equations (2.2) and (2.5) with Equations (3.9) and (3.14), respectively, the computational procedures of MQCB developed in the literature can be used to simulate molecular dynamics including vibration, rotation and spin motions. Before this new QM/MM approach becomes workable, we have to verify that the governing equations derived in Section 3 yield correct molecular dynamics. The verification is based on the comparison with the quantum mechanical description of a diatomic molecule for which the Schrödinger equation relating to the nuclear motion has an analytical solution. The test model is illustrated in Figure 1a, which is an equivalent single-particle model of a diatomic molecule with rotational and vibrational motion described by the spherical coordinates. The effective internuclear potential is modeled by the Morse function and the corresponding Schrödinger equation is solved analytical in the appendix with eigenfunction Ψ , , ( , , ) = , ( )Θ , ( )Φ ( ) given by Equation (A15) and eigenvalue , given by Equation (A17). All of the parameters and constants appearing below refer to the Appendix.
As a comparison, if the standard Bohmian mechanics based on Equation (2.2) is used to simulate the molecular dynamics in the eigenstate Ψ , , ( , , ), we will find the molecule to be motionless in all the states with = 0. The constitution of the complex momentum in Equation (3.5) explains why it is possible to observe non-vanishing momentum in cases where the Bohmian momentum ∇ vanishes.
The first test on the accuracy of the state-dependent molecular dynamics is to examine the quantization laws existing in the eigenstate Ψ , , ( , , ). In the previous section, we have seen that in the conventional MM description, the three conservation constants depend continuously on the initial conditions, being unable to manifest the quantization laws. Now replacing Equations (2.8) and (2.9) by Equations (3.8) and (3.9), we can derive the expected quantization laws. With the eigenfunction Ψ , , ( , , ) given by Equation (A15) and ( , , ) given by Equation (3.10), the Hamiltonian , the squared angular momentum and the z-component angular momentum defined in Equation (3.8) can be computed as: We find that the resulting values of , and are independent of time and are quantized according to the three quantum numbers , and . Unlike the probabilistic interpretation in standard quantum mechanics, here, we have given a dynamic interpretation of the quantization laws by showing that the three physical quantities have constant discrete values along any quantum trajectory ( ( ), ( ), ( )) solved from Equation (4.1) in the quantum state specified by the three quantum numbers ( , , ).
The second test is on the consistence between the predictions of the equilibrium bond length made between quantum mechanics and state-dependent molecular dynamics. We will show that the equilibrium bond length that maximizes the radial probability satisfies the dynamical equilibrium condition / = 0. We first prove this property for the ground vibrational state. Referring to the Appendix, the corresponding radial wave function with = 0 is given by . (4.4) The above radial dynamics has a stable equilibrium point at = , i.e., (4.5) The other equilibrium point is at = 0, i.e., ̅ = ∞, corresponding to the condition of molecular dissociation. Equation (4.5) expresses the equilibrium bond length ̅ as an explicit function of the angular momentum quantum number J. This closed-form expression describes analytically how the equilibrium bond length ̅ increases monotonically with the angular quantum number J. Of significance is that the equilibrium bond length ̅ obtained from the molecular dynamics (4.1a) always coincides with that obtained by QM method. This coincidence has its theoretic origin. We recall the definition of the radial probability: , (4.6) by noting that , ( ̅ ) given by Equation (A13) is a real function of ̅ . On the other hand, according to the dynamic Equation (4.1a), the equilibrium position ̅ satisfies the condition: . (4.7) It appears that the equilibrium condition (4.7) is just the condition requiring the radial probability , in Equation (4.6) to have an extreme value. All of the properties obtained from the dynamics equation (4.1a) can be re-derived from the action of the radial total force ̅ ( ̅ ). As an illustration, we consider the case with = 0 for which Θ , ( ) = 1 and the total potential defined in Equation (3.12) has a simple expression, , where we note = 2 − 1 and = from Equation (A14) for the case of = 0. The radial total force ̅ ( ̅ ) now can be determined from ( ̅ ) as . (4.9) The internuclear distance free from radial force can be found from the condition ̅ ( ̅ ) = 0, which leads to . (4.10) This value is exactly the equilibrium bond length ̅ already obtained in Equation (4.5) by using the equilibrium condition ̅ / = 0. The evaluation of ̅ ( ̅ ) in the vicinity of ̅ gives ̅ ( ̅ ) < 0, if ̅ > ̅ , and ̅ ( ̅ ) > 0, if ̅ < ̅ . The sign of ̅ ( ̅ ) in the neighborhood of ̅ indicates that the two nuclei attract each other when their distance is longer than ̅ and repel each other when their distance is shorter than ̅ . Consequently, there is always a restoring force to make ̅ return to its equilibrium position ̅ .

Equilibrium positions
The complex trajectories of ̅ ( ) solved from Equation (4.4) for the four types of diatomic molecules are shown in Figure 3a. It can be seen that the complex trajectories are closed contours circulating about their respective equilibrium positions ̅ computed from Equation (4.5) by using the molecular parameters listed in Table 1. We find that the period of vibration is independent of the actual trajectories and is quantized with respective to the angular quantum number J. This trajectory-independent property can be proven by applying the residue theorem to Equation (4.4): , (4.11) where the contour is an arbitrary closed trajectory solved from Equation The values of the four molecules and their periods of vibration computed by Equation (4.12) are listed in Table 1. Also shown in Table 1 are the experimental values of the ground-state vibration frequencies. The comparison between the computational and experimental results gives a relative error about 1%, which is caused by the inexistence of an exact solution to the Schrödinger Equation (A1). The radial wave function , given by Equation (A13) is only a second-order approximation as described by Equation (A6). The approximation gets better as is closer to the equilibrium bond length of the Morse potential. The vibration period given by Equation (4.12) can also be derived by using a state-dependent force constant. A quantum force constant is defined as the second-order derivative of the total potential ( ̅ ) evaluated at the equilibrium position ̅ : . (4.13) The resulting force constant is state-dependent by noting that depends on the three quantum numbers ( , , ) as given by Equation (3.12). A specific has been expressed explicitly in Equation (4.8) for the case of = = 0, and the substitution of into Equation (4.13) yields the force constant as: . (4.14) Then the classical relation between the vibration period and the force constant gives: , (4.15) which reproduces the result of Equation (4.12) derived by state-dependent molecular dynamics. The above discussions of ground vibrational state reveal a remarkable observation that by replacing the internuclear potential with the total potential , the conventional MM turns out to be state-dependent MM, which then yields consistent results with QM. This observation applies to every quantum state of the molecule. In every quantum state Ψ , , , there is a total potential , which completely determines the molecular quantum dynamics within this state. Figure 4 plots the total potentials for the first three vibrational states and their accompanying rotational states. It can be seen that for a vibrational state with quantum number , there are + 1 sub-shells, within which rotational states with different quantum number J reside. In excited vibrational states ≥ 1, the quantum dynamics (4.1a) has multiple equilibrium points, and the total potential possesses multiple-shell structure. To illustrate this property, we consider the first excited state as an example. The corresponding radial wavefunction is given by: , (4.16) and the radial dynamics is obtained from Equation (4.1a) as: . (4.17) There are two stable equilibrium points in this state, which can be solved from the condition / = 0 as: ( ) On the other hand, the total potential in the first excited state is given by Equation (3.12): . (4.19) By evaluating the constant at = 1 with = 2 and = 3, we have: . (4.20) Substituting this into Equation (4.19), the variation of ( ̅ , ) with respect to the angular quantum number J is plotted in Figure 4. It appears that there is a node at = + 1, where the total potential approaches to infinity. This infinite potential barrier separates ( ̅ , ) into two shells. The lowermost point of the inner shell corresponds to the inner equilibrium point ̅ ( ) , while the lowermost point of the outer shell corresponds to the outer equilibrium point ̅ ( ) , as given by Equation (4.18). Also shown in Figure 4 is the total potential of the second excited state, which has three shells separated by the two nodes.

State-Dependent Orbital and Spin Dynamics
So far our discussions on molecular quantum dynamics are restricted to the radial vibration motion governed by Equation (4.1a). To take angular dynamics into account, all the three equations in Equation (4.1) have to be considered. In the previous sections we have seen that the QM descriptions of the quantum state Ψ , , can be reproduced in terms of the state-dependent molecular dynamics. Here, we go one-step further to show that the spin dynamics, which cannot be described by the spatial wave function Ψ , , , emerges naturally from the molecular dynamics where and are given by Equation (A14) with = 0. We can see that even though the orbital angular quantum number J is set to zero, the angular velocity is actually not zero. The main reason why spin motion cannot be attributed to orbital motion in standard quantum mechanics is due to the definition of the orbital angular momentum . To make this point more apparent, we inspect the expression of defined in Equation (1) (2) eq eq 2 3 8 9 , 2 2 z z α α In the real domain, the only solution is = 0 and = /2, which is the well-known spinless motion. However, in the complex domain, Equation (5.3) has a nontrivial solution: which is just the dynamics in Equation (5.1) in dimensionless form. In other words, the condition of = 0 does not completely nullify the complex-valued angular motion, and we can regard the spin dynamics as the remnant angular dynamics as the orbital angular momentum is zero.
The integration of the dynamics with = + gives the three regions of trajectories as shown in Figure 5a. Within the Ω region, the sign of / changes alternatively and produces zero mean velocity. Within the Ω region, ( ) is monotonically decreasing with a mean value of / equal to −1. Within Ω region, ( ) is monotonically increasing with a mean value of / equal to one. Although orbital angular momentum is completely depressed by the condition = = 0, non-zero angular motion does exist in the regions of Ω and Ω . To relate this remnant angular motion to spin, the next step is to identify the magnitude of this remnant angular momentum. In the regions of Ω and Ω , we note the following relations for the complex variable = + : , Therefore, angular momentum in the direction becomes: In the region Ω , the ( ) dynamics exhibits periodic motion around the equilibrium points and yields zero mean angular velocity. According to the behavior of the dynamics, three spin regions can be defined, The spin dynamics we discuss so far originates from the wavefunction Ψ = , ( )Θ , ( )Φ ( ) with Θ , ( ) given by (cos ), the first-type Legendre function. Indeed, (cos ) is only a special solution for Θ , ( ) whose general solution can be represented by: where ( ) is the second-type Legendre function. It can be shown by integration Equation (4.1) with Θ , ( ) = (cos ) that the direction of the angular motion produced by (cos ) is always anti-parallel to that produced by (cos ). Therefore, in the same quantum state specified by the spatial quantum numbers ( , , ), the molecule can behave either according to the spin dynamics with Θ , ( ) = (cos ) or to the anti-spin dynamics with Θ , ( ) = (cos ). For the former case, we say that the molecule is in a sub-state of ( , , ) denoted by ( , , , 1/2), and for the latter case, the molecule is in the sub-state ( , , , −1/2). Furthermore, an entangled spin dynamics can be simulated according to the state-dependent molecular dynamics by expressing Θ , ( ) as a linear combination of (cos ) and (cos ) as described by Equation (5.8).
In the ground state, we have seen that the orbital angular motion vanishes and the remnant angular motion in the direction emerges as the spin dynamics. We proceed to demonstrate that in excited states, orbital and spin motions coexist and both contribute to angular momentum. To examine the interaction between orbital and spin dynamics, we next consider quantum states with = 1 and = 0. Since the -component angular momentum is zero in these quantum states, we are interested in where the orbital angular momentum emerges and how it interacts with the spin angular momentum. The and dynamics for these quantum states are governed by Equations (4.1b) and (4.1c) as: . (5.9) The trajectories of ( ) on the complex − plane are illustrated in Figure 5b, where it can be seen that as in the ground state, three regions representing three angular motions come about. In the region Ω ( ≫ 0), is monotonically decreasing, while in the region Ω ( ≪ 0), is monotonically increasing. In the central region Ω , is a periodic function of time. In conjunction with the relation (5.5), the mean angular momentum of Equation (5.9) becomes: )   2   cos  sin  ,  0,  2  sin  cos  2  cos  cos  sin  2 sin  ,  0,  2  sin  cos  2   I   I   I   I   R  R  I Comparing with Equation (5.7), we find that in the regions of Ω and Ω , the angular momentum contains an additional component ℏ contributed from the orbital motion = 1 apart from the spin angular momentum ℏ/2, indicating that the orbital angular motion, in the states with = 1, and = 0, is produced solely by the dynamics because of = 0. For a quantum state with quantum number ≠ 0 and ≠ 0 , the and dynamics can be expressed by: Making use of Equation (5.5) once again, we obtain a general expression for the angular momentum in the and directions as: In summary, in the Ω region, one has zero-mean dynamics and a quantized -component angular momentum ℏ, which is a well-known result in QM; while in the Ω and Ω regions, one has an angular momentum ±( + 1/2)ℏ in the direction and a zero-mean orbital angular momentum in the direction, which is otherwise unknown to QM. The angular momentum given by Equation (5.12), which is derived from the state-dependent molecular dynamics, provides us with a trajectory-based method to determine the rotational energy of a diatomic molecule and its rotational spectrum. The most used formula of rotational energy takes the following form: , (5.13) which treats the molecule as a rigid rotor with moment of inertia = /2. A typical rotational spectrum consists of a series of peaks corresponding to transitions between adjacent levels satisfying ∆ = ±1: 14) The rigid-rotor model assumes a constant internuclear distance and neglects the stretch of the bond as the molecule rotates. To account for bond stretching due to rotation, we consider a refined rotation-energy component in the eigenvalue expansion given by Equation (A20), from which the wavenumber of the transition from + 1 to can be derived as: where is a correction factor known as the centrifugal distortion constant. Table 2 lists the measured far infrared absorption spectrum ( ) [29] for the diatomic molecule and the predicted spectrum by the rigid model ( ) and the non-rigid model ( ) with = 10.58 cm and = 5.6 × 10 cm determined from the molecular parameters listed in Table 1. Instead of using the quantum rotational energy obtained from the eigenvalues of Schrödinger equation, we will give an estimate of the wavenumber of transition based on the classical expression of the rotational energy = with and computed by the state-dependent molecular dynamics: , (5.16) where is given by Equation (5.12) and the bond length is denoted by to emphasize its dependence on the angular quantum number J as given by Equation (4.5) and Equation (A14). An alternative estimate of the wavenumber caused by rotation absorption transition turns out to be: . (5.17) With the molecular parameters of given by Table 1, the bond length is calculated by using Equation (4.5) and then substituted into Equation (5.17) to yield the estimated spectrum . Table 2 compares the four predictions  ,  , and with the measured rotational spectrum . It can be seen that the wavenumber is closest to , because is obtained by tuning the parameters and , so that Equation (5.15) best fits the experimental data. The curve fitting gives the best values as * = 10.44 cm and * = 5.2 × 10 cm . The two wavenumbers and can be compared on the same footing, since both adopt the same molecular parameters listed in Table 1 in the computation. Their difference lies in the model of rotational energy used to compute wavenumbers. The wavenumber is based on the eigen-energy model (A20), while is based on the QHM model (5.16). Theoretically, the description of state-dependent molecular dynamics by QHM is equivalent to QM; however, the actual accuracy depends on the degree of approximation involved in the solution to the wave function. The equilibrium bond length computed from Equation (4.5) is based on the radial wave function , given by Equation (A13), which is an approximate solution to the radial Schrödinger Equation (A5) by employing the second-order expansion (A6). When the quantum number J increases, the bond length deviates further from the equilibrium position of the Morse potential, and the accuracy of the expansion (A6) gets worse. This tendency explains the degradation of the prediction accuracy of with respect to the measured rotational spectrum .

Conclusions
A new QM/MM approach called quantum Hamilton mechanics (QHM), is proposed in this paper to establish state-dependent molecular dynamics (SDMD) in such a way that the governing equations of SDMD can be derived by MM with solutions compatible with QM. As a complex extension of Bohmian mechanics, QHM is coordinate-independent and especially suitable in curvilinear coordinates to simulate coupled orbital/spin dynamics. The correctness of SDMD has been verified by comparing with the quantum mechanical description of a diatomic molecule for which Schrödinger equation has an analytical solution. The resulting SDMD simultaneously satisfies the continuous-time dynamics governed by MM and the quantized dynamics governed by QM. QHM can be incorporated into the current framework of the mixing quantum/classical Bohmian (MQCB) method to simulate molecular dynamics. The incorporation of QHM with MQCB enables a trajectory interpretation of orbital-spin interaction and makes it possible for us to simulate spin entanglement in molecular dynamics.