# DECIDE: A Deterministic Mixed Quantum-Classical Dynamics Approach
## 4 Hydrogen Bond Model
### Model  
The Weyl-ordered, partially Wigner-transformed Hamiltonian of the hydrogen bond model is given by  
\[\hat{H}_{W}=\frac{\hat{p}^{2}}{2m}+\frac{P^{2}}{2M}+A^{0}\hat{x}^{2}+\frac{1}{2} k_{1}X\hat{x}+\frac{1}{2}k_{1}\hat{x}X, \tag{19}\]  
where \(\hat{x}/\hat{p}\) is the position/momentum operator of the proton, \(X/P\) is the position/momentum of the solvent coordinate, and \(A^{0}\) and \(k_{1}\) are constants. Plots of the protonic potential for two values of \(X\) are shown in Figure 2.  
A convenient basis for solving this proton transfer problem is a set of quantum harmonic oscillator wave functions, i.e.,  
\[\phi_{n}(x)=\langle x|\phi_{n}\rangle=(2^{n}n!\sqrt{\pi})^{-1/2}b^{1/2}H_{n}(bx )\exp[-b^{2}x^{2}/2], \tag{20}\]  
where \(H_{n}\) is the \(n\)th Hermite polynomial, \(b\) is an arbitrary constant, and \(n=0,1,2,\ldots\). Note that if \(X=0\) in Equation (19), then the exact solution of the time-independent Schrodinger equation for this model is \(\phi_{n}(x)\) with \(b=\left(2 mA^{0}/\hbar^{2}\right)^{1/4}\). In our simulations, we take \(b=\left(8 mA^{0}/\hbar^{2}\right)^{1/4}\). The eigenfunctions of \(\hat{H}_{W}\) may be expanded in this basis as \(|\psi_{i}\rangle=\sum_{i}c_{n}^{i}|\phi_{n}\rangle\), where \(c_{n}^{i}\) is the coefficient of the \(n\)th basis function. For this expansion, we use 12 basis functions (\(n=0,\ldots,11\)). Substituting the eigenfunction expansion into the time-independent Schrodinger equation, we obtain an eigenvalue problem of the form \(Hc=cE\), where \(H\) is the Hamiltonian matrix expressed in the chosen basis, \(c\) is a matrix containing the \(c_{n}^{i}\)'s of the \(i\)th eigenfunction, and \(E\) is a diagonal matrix containing the eigenvalues \(E_{i}(X)\). The Hamiltonian matrix elements are calculated by numerical integration, and the eigenvalue problem is solved (to obtain the \(c_{n}^{i}\)'s and eigenvalues) using built-in functions in MATLAB. Figure 3 shows the ground- and first-excited state adiabatic potential energy surfaces [\(E_{0}(X)\) and \(E_{1}(X)\), respectively] calculated using this method.  
It turns out that the time-independent Schrodinger equation for the Hamiltonian in Equation (19) can be solved analytically (and exactly) by first rewriting the potential energy as \(A^{0}(x+\frac{k_{1}}{2A^{0}}Xx)^{2}-\frac{k_{1}^{2}}{4A^{0}}X^{2}\). The resulting adiabatic energies are  
\[E_{n}(X)=\sqrt{\frac{2A^{0}}{m}}(n+\frac{1}{2})\hbar-\frac{k_{1}^{2}}{4A^{0}} X^{2},\ \ \ n=0,1,2,3,\ldots \tag{21}\]  
Figure 2: Plots of the protonic potential along the proton coordinate \(x\). The blue and red curves represent the potential for two different values of the solvent coordinate \(X\), which modulates the position and depth of the protonic potential. The dashed line corresponds to the ground state energy for a particular value of the solvent coordinate.  
where \(n\) labels the energy level. Comparing our variational results to the exact ones, we find that good agreement is obtained when \(X\) is close to 0, but differences emerge when \(|X|>0.25\) (see Figure 3). The convergence of the variational result can be improved by increasing the number of basis functions; however, the level of accuracy achieved with 12 basis functions is sufficient for our purposes.