% REVTeX 3.0
% For APS galley-mode format:
%\documentstyle[pra,aps]{revtex}
\documentstyle[pra,aps,twocolumn]{revtex}
%\documentstyle[eqsecnum,aps,amsfonts]{revtex}
%
% For preprint format:
%\documentstyle[preprint,pra,aps]{revtex}
%\documentstyle[preprint,eqsecnum,aps,amsfonts]{revtex}
\tighten
%
\begin{document}
% \draft command makes pacs numbers print
\draft
\title{Self-consistent-field perturbation theory of
molecular vibrations}
% repeat the \author\address pair as needed
\author{Alexei V. Sergeev\cite{vavilov} and David Z. Goodson}
\address{Department of Chemistry, Southern Methodist University,
Dallas, Texas 75275}
\date{\today}
\maketitle
\begin{abstract}
Perturbation theory is used to perform noniterative
calculations of energy eigenvalues of the coupled ordinary
differential equations that result from imposing
separability assumptions in terms of normal coordinates
on vibrational wavefunctions. Various model Hamiltonians
with 2 or 3 coupled normal modes are studied and the
increase of computational cost with the number of
degrees of freedom is analyzed. Quadratic
Pad\'e approximants of the perturbation expansions are
rapidly convergent, and directly yield complex numbers for
resonance eigenvalues. For a 3-mode system, results are
obtained within partial separability assumptions, with
a pair of modes left coupled. Large-order perturbation
theory with partial separability is suggested as an
alternative to low-order exact perturbation theory.
\\
\end{abstract}
% insert suggested PACS numbers in braces on next line
\pacs{}
% body of paper here
\narrowtext
\section{Introduction}
The vibrational self-consistent-field (VSCF) approximation
\cite{bowmanrev,carneyrev,gerberrev},
in which the molecular vibrational wavefunction is expressed as
a separable product of normal-coordinate functions, extends
normal-mode analysis to mildly coupled systems with only
a modest increase in computational difficulty. For tightly
bound molecules such as H$_2$O and CO$_2$ the separability
assumption introduces a systematic error in the energy on the
order of 1\%. However, it makes the analysis of
large systems computationally feasible \cite{bpti,argon,jelski}.
To improve the accuracy, VSCF wavefunctions can be used
as the basis for a ``state-interaction'' calculation,
analogous to the configuration-interaction method for
electronic structure calculations, but this greatly increases
the computational cost.
Here we study VSCF approximations in the context
of large-order perturbation theory.
It was recently demonstrated that large-order
Rayleigh-Schr\"odinger perturbation theory
can have a significant advantage in computational cost over
variational methods for solving the
Schr\"odinger equation for vibration spectra \cite{anharm}
and vibration-rotation spectra \cite{roviblet}. In addition
it is a convenient approach for calculating
resonance energies, because the computational algorithm for
resonance eigenvalues is equivalent to that for bound states,
as long as appropriate summation approximants are used
\cite{popovstark,germann,sergeev}.
In Ref.~\cite{gscf} a general algorithm was developed
for solving SCF equations with perturbation theory.
Perturbation theory decouples the equations order by order,
thus yielding direct, noniterative, solutions for asymptotic
expansions of wavefunctions and energies. Furthermore, it was
shown that perturbation theory was straightforward to apply with
partial separability, with the wavefunction expressed
as a separable product of functions that depend on clusters of
normal coordinates.
We have several motivations for developing an SCF perturbation
theory of molecular vibrations. The first is to provide a
feasible strategy for applying perturbation theory to large
molecules. The cost of computing high-order
coefficients of perturbation expansions increases rapidly
with the number of degrees of freedom. Large-order SCF
perturbation theory is an alternative to low-order exact
perturbation theory. Thus one can eliminate
the error from premature truncation of the
expansion at the cost of introducing a systematic error
from a separability assumption. Another motivation
is to study the accuracy of partial separability assumptions.
Partial separability offers an alternative to the
state-interaction method for systematically improving the
accuracy of VSCF calculations. It can also be used as
a probe of the degree of interaction between particular
normal modes. Finally, SCF perturbation theory can be applied
to a wide variety of problems including, for example, the
electronic Schr\"odinger equation \cite{gscf,dzgjcp,scfmag}.
The molecular vibration problem is a straightforward test case
for the general theory.
In Section~\ref{sec:method} we outline the calculation method.
In Section~\ref{sec:results} we briefly treat a special case
in which the SCF equations can be solved exactly, and then
present results from perturbation theory for various 2- and
3-mode systems that have been treated previously in the
literature. The emphasis is on determining the rate of
convergence of the perturbation expansions and on studying
the accuracy of the SCF approximation for imaginary parts of
resonance energies. For a 3-mode system we compare the accuracy
of partial separability approximations to that of
state-interaction calculations. Section~\ref{sec:cost} contains
an analysis of how the computational cost of large-order
SCF perturbation theory increases with the number of degrees
of freedom. In Section~\ref{sec:discussion} we discuss potential
advantages of this theory and suggest future applications.
\section{SCF theory}
\label{sec:method}
Consider a Hamiltonian of the form
\begin{equation}
H=w({\bf q})+\sum_{\alpha=1}^N h^{(\alpha)}(q_\alpha),
\label{genlham}
\end{equation}
where
\begin{equation}
h^{(\alpha)}=-{1\over 2}{\partial^2\over\partial q^2_\alpha}
+{1\over 2}\omega_\alpha^2q^2_\alpha
+\eta_{\alpha,1}\, q^3_\alpha+\eta_{\alpha,2}\, q^4_\alpha+\cdots.
\label{halpha}
\end{equation}
${\bf q}=\{q_1,q_2,\dots,q_N\}$ is the set of harmonic normal
coordinates, which diagonalizes the harmonic expansion of $H$,
and $w({\bf q})$ is a polynomial with degree of at least 3
that couples the $q_\alpha$.
This type of Hamiltonian could arise from the use of a
finite-degree polynomial fit to {\it ab initio} or empirical
information for the potential energy operator or from the use
of a Taylor expansion of an analytic potential function.
The SCF method consists of
replacing the partial differential Schr\"odinger equation
$H\Psi=E\Psi$ with a set of coupled ordinary differential equations,
\begin{mathletters}
\begin{equation}
\left( H^{(\alpha)}
-\epsilon^{(\alpha)}\right)\psi^{(\alpha)}(q_\alpha)=0,
\label{diffeqs}
\end{equation}
\begin{equation}
H^{(\alpha)}=
h^{(\alpha)}(q_\alpha)+\bar{w}^{(\alpha)}(q_\alpha),
\label{Halpha}
\end{equation}
\begin{equation}
\bar{w}^{(\alpha)}(q_\alpha)=
\bigg\langle \prod_{\beta\ne\alpha}\psi^{(\beta)}(q_\beta)\bigg|
w({\bf q})\bigg|
\prod_{\beta'\ne\alpha}\psi^{(\beta')}(q_{\beta'})\bigg\rangle.
\label{wbar}
\end{equation}
\label{vscfeqs}
\end{mathletters}
\noindent
Equations\ (\ref{vscfeqs}) result from substituting for the
wavefunction a separable product
$\Psi=\prod_{\alpha}\psi^{(\alpha)}(q_\alpha)$ and then applying
the variational principle for the energy with $\Psi$
constrained to be normalized.
Let
\begin{equation}
H=
\sum_{\alpha=1}^N\left(
-{1\over 2}{\partial^2\over\partial q_\alpha^2}
+{1\over 2}\omega_\alpha^2q_\alpha^2\right)
+H'
\label{hamiltonian}
\end{equation}
where $H'$ is the sum of $w$ and all
of the anharmonic terms in the $h^{(\alpha)}$.
We will write $H'$ in the form
\begin{equation}
H'=\sum_{k=1}\kappa^{k/2}\left[
w_k({\bf q})+\sum_\alpha \eta_{\alpha,k}q_\alpha^{k+2}
\right] ,
\label{Hprime}
\end{equation}
where $\kappa$ is a perturbation parameter and the $w_k$ are
sums of monomials of degree $k+2$ that couple normal
coordinates. The physical solution corresponds to $\kappa=1$.
It follows that $\psi^{(\alpha)}$ and the $\epsilon^{(\alpha)}$
have asymptotic expansions
\begin{equation}
\psi^{(\alpha)}(q_\alpha)
=\sum_{k=0}^\infty\kappa^{k/2}\psi_k^{(\alpha)}(q_\alpha),
\label{psix}
\end{equation}
\begin{equation}
\epsilon^{(\alpha)}=
\sum_{k=0}^\infty\kappa^{k/2}\epsilon_k^{(\alpha)}.
\label{epsilonx}
\end{equation}
Substituting these into Eqs.~(\ref{vscfeqs}) and
collecting terms multiplied by $\kappa^{k/2}$ yields a set of
$N$ equations for each value of $k$. These
equations are coupled implicitly through the dependence of
$\bar{w}^{(\alpha)}$ on the $\psi^{(\beta\ne\alpha)}$.
The total vibrational energy is
\begin{eqnarray}
E&=&{1\over N}\sum_\alpha\epsilon^{(\alpha)}
+{N-1\over N}\sum_\alpha\left\langle h^{(\alpha)}\right\rangle
\nonumber\\
&&\quad
+\;{N-2\over N}\sum_{\alpha_1<\alpha_2}
\left\langle w^{(\alpha_1,\alpha_2)}\right\rangle
+\;\cdots
\nonumber\\
&&\quad
+\;{1\over N}\sum_{\alpha_1<\alpha_2<\cdots<\alpha_{N-1}}
\left\langle w^{(\alpha_1,\alpha_2,\dots,\alpha_{N-1})}\right\rangle,
\label{totalE}
\end{eqnarray}
\noindent
where the $w^{(\alpha_1,\dots)}$ are the terms in $w({\bf q})$
that couple the designated normal coordinates.
Here we will only consider functions $H'$ such that
$H'(-{\bf q},\kappa^{1/2})=H'({\bf q},-\kappa^{1/2})$.
It follows that the asymptotic expansion of $E$ will only
include integer powers of $\kappa$.
Now consider the expansions in powers of $\kappa^{1/2}$ that
result for the terms
$\bar{w}^{(\alpha)}(q_\alpha)$. Substituting Eq.~(\ref{psix})
into Eq.~(\ref{wbar}) yields a sum of terms
$\big\langle \psi_{k_1}^{(\beta)}\big| w_{k_2}
\big| \psi_{k_3}^{(\beta)} \big\rangle$. For terms at order
$\kappa^{k/2}$, we have $k_1+k_2+k_3=k$. However, the
smallest possible value of $k_2$ is 1, which implies that the
largest possible value for $k_1$ or $k_3$ is $k-1$.
Therefore, the equation that determines $\psi_k^{(\alpha)}$
depends on the $\psi_{k'}^{(\beta\ne\alpha)}$ only for
$k' {\sqrt{3}\over 9}
\left(n_y+{1\over 2}\right)^{-1}\omega_x^2\omega_y^3 ,
\label{omegacondition}
\end{equation}
in which case the eigenstate represents a resonance.
The energy of the resonance is $\Re E$ and the
width is $\Gamma\approx 2|\Im E|$ \cite{connor}.
Since the value of the last term in Eq.~(\ref{Ecase1})
does not depend on $n_x$, the width of the resonance,
within the SCF approximation, depends only on $n_y$.
Figure~\ref{fig1} compares the SCF results to exact
results for various eigenstates. (The ``exact'' energies
were calculated with large-order perturbation theory
\cite{anharm} without the separability assumption.) We have
subtracted $(n_x+{1\over 2})\omega_x$ from the real part
in Fig.~\ref{fig1} to make the SCF curves independent
of $n_x$. The exact results show a weak dependence
on $n_x$. The exact imaginary part is nonzero for all
$\lambda>0$, although it becomes exponentially small as
$\lambda$ approaches zero.
Note that the accuracy of the SCF approximation for the state
$(n_x,n_y)=(0,1)$ is greater than that for (1,0). This is
perhaps due to the fact that Eq.~(\ref{Ecase1}) includes
nonlinear behavior only in $n_y$; the $n_x$ spectrum within
the SCF approximation is purely
harmonic. We find that the SCF results for the $(0,n)$
states in general tend to be more accurate than those for
the $(n,0)$ states, although this trend becomes less noticeable
for very large $n$.
\subsection{2-mode systems}
\label{sec:2mode}
Consider the Hamiltonian
\begin{eqnarray}
H&=&-{1\over 2}{\partial^2\over\partial x^2}
-{1\over 2}{\partial^2\over\partial y^2}
\nonumber\\
&&\quad
+{1\over 2}\omega_x^2x^2\! +{1\over 2}\omega_y^2y^2
-\lambda (xy^2+\eta x^3),
\label{Hcase2}
\end{eqnarray}
\noindent
which was studied by Eastes and Marcus \cite{eastes} for
various values of $\omega_x$, $\omega_y$, $\lambda$ and $\eta$.
Let $\kappa^{1/2}=\lambda$ be the perturbation parameter.
Expansions for the exact energy eigenvalues,
without separability, were calculated in Ref.~\cite{anharm}.
Although these expansions are in principle divergent if
taken to high enough order \cite{cizek}, in practice the partial
sums were rapidly convergent and the rate of convergence could be
increased by Pad\'e summation. We find that the SCF expansions
behave similarly, but the rate of convergence is even more rapid
than for the exact perturbation expansions. Furthermore, we find
that the convergence for both the exact and the SCF expansions
can be improved by using quadratic Pad\'e approximants
\cite{shafer},
\begin{equation}
S_{[n,m,l]}={1\over 2R_n}
\left( Q_m\pm\sqrt{Q_m^2+4P_lR_n}\right),
\label{quad}
\end{equation}
where $P_l$, $Q_m$, and $R_n$ are polynomial of degrees
$l$, $m$, $n$, respectively, in the expansion parameter.
These polynomials are determined from the asymptotic relation
\begin{equation}
RE^2+QE-P\sim O\big(\kappa^{n+m+l+2}\big),
\label{quadrelation}
\end{equation}
with $E$ in Eq.\ (\ref{quadrelation}) representing the
asymptotic expansion of the energy. A comparison of the
convergence of the quadratic approximants for the exact and
SCF expansions is given in Fig.~\ref{fig2}.
Quadratic approximants are especially useful for resonance
eigenvalues \cite{popovstark,germann,sergeev}. Since the expansion
coefficients that we calculate for $E$ are real numbers,
the standard linear
Pad\'e approximants $S_{[n,m]}=P_n/Q_m$, defined by
$QE-P\sim O(\kappa^{n+m+1})$,
will always be real. However, quadratic approximants will
give an imaginary part when the discriminant is negative.
Figure~\ref{fig3} compares the converged values of SCF and
exact energies. The qualitative trends in the accuracy of the
SCF approximation are similar to those in Fig.~\ref{fig1}.
Another model for vibrational resonances is
\begin{mathletters}
\begin{equation}
H=-{1\over 8}{\partial^2\over\partial x^2}
-{1\over 8}{\partial^2\over\partial y^2}+V(x,y),
\label{Hcase3}
\end{equation}
\begin{equation}
V(x,y)=
{1\over 2}\omega_x^2x^2e^{-x^2}+{1\over 2}\omega_yy^2
- {1\over 2}\omega^2_yx^2y^2e^{-x^2},
\label{Vcase3}
\end{equation}
\label{case3}
\end{mathletters}
\noindent
which was devised by Waite and Miller \cite{waite} to describe
unimolecular dissociation, and has
been studied within the SCF approximation
by Christoffel and Bowman \cite{christoffel}
and by Farrelly and coworkers \cite{reinhardt,farrelly}.
The potential has a symmetric double barrier in the $x$ mode.
It is convenient to express this problem as a
semiclassical perturbation theory \cite{mlodinow}, with the
perturbation parameter multiplying the kinetic energy:
\begin{equation}
\left[\kappa^2\left( -{1\over 8}{\partial^2\over\partial x^2}
-{1\over 8}{\partial^2\over\partial y^2}\right)
+V(x,y)-E\right]\Psi=0.
\label{HPTcase3}
\end{equation}
If we make the substitutions $x\to\kappa^{1/2}x$,
$y\to\kappa^{1/2}y$, and $E\to\kappa E$ and then expand
the potential in powers of $\kappa^{1/2}$, then we obtain
a Hamiltonian that is essentially in the form given by
Eqs.~(\ref{hamiltonian}) and (\ref{Hprime}).
SCF results from quadratic Pad\'e summation are presented in
Table~\ref{tab1}. Figure~\ref{fig4} shows representative
convergence plots. The accuracy obtained at given order for the
imaginary part is approximately the same as that for the real part.
\subsection{A 3-mode system}
\label{sec:3mode}
Consider the Hamiltonian
\begin{eqnarray}
H&=&-{1\over 2}{\partial^2\over \partial x^2}
-{1\over 2}{\partial^2\over \partial y^2}
-{1\over 2}{\partial^2\over \partial z^2}
\nonumber\\
&&+{1\over 2}\omega_x^2x^2
+{1\over 2}\omega_y^2y^2
+{1\over 2}\omega_z^2z^2
\nonumber\\
&&-\kappa^{1/2}(xy^2+yz^2+\eta x^3 +\eta y^3),
\label{Hcase4}
\end{eqnarray}
\noindent
with $\omega_x=0.7$, $\omega_y=1.3$, $\omega_z=1.0$,
$\eta=0.1$, and $\kappa^{1/2}=0.1$. This is another standard
test case \cite{noid,bowmancpl,moelples}. There are now
four possible normal-mode separability assumptions ---three
in which a pair of modes is left coupled and one
with all modes separated.
Christoffel and Bowman \cite{bowmancpl} used an iterative
method to calculate SCF energies for this system and then
used an expansion of the exact wavefunction in terms of a finite
set of SCF wavefunctions to obtain a variational
estimate for the exact energy. This ``state-interaction''
method will converge to the exact energy (for bound states)
as the size of the SCF basis is increased. The computational
cost is not too high for 3-mode systems, but increases rapidly
with the number of degrees of freedom. Norris {\it et al.}
recently studied this system using M\o ller-Plesset
perturbation theory \cite{moelples}.
We consider the
cases $\Psi_{\{x;y,z\}}=\psi^{(a)}(x)\psi^{(b)}(y,z)$ and
$\Psi_{\{z;x,y\}}=\psi^{(a)}(x,y)\psi^{(b)}(z)$.
(The case $\{y;x,z\}$ is for this system
equivalent to full separability, since none of the coupling
terms in Eq.~(\ref{Hcase4}) couple $x$ and $z$.)
Figure~\ref{fig5} compares the accuracy of the various
separability assumptions for the states ($n_x$,0,0) and
(0,$n_y$,0). As $n_x$ increases, the accuracy of $\{z;x,y\}$
holds steady while $\{x;y,z\}$ and $\{x;y;z\}$
become less accurate. Increasing $n_x$ causes
an increase in the expectation value of $|x|$, which increases
the importance of the term $xy^2$ in Eq.~(\ref{Hcase4});
hence, the approximations with $x$ and $y$ separated
become less accurate as this term becomes more
important. As $n_y$ increases, the $\{x;y,z\}$ separability
approximation becomes the most accurate. This appears to be
analogous to the trend seen in Fig.~\ref{fig1} for the
exactly solvable problem: If the coupling between two normal
coordinates has the form $q_1q_2^2$, then the SCF
approximation tends to be more accurate for $(0,n_2)$
excitations than for $(n_1,0)$ excitations. Thus,
$\{x;y,z\}$ separability with $xy^2$ coupling and $n_y$
excitation tends to be more accurate than $\{z;x,y\}$
with $yz^2$ coupling and $n_y$ excitation.
Also shown in Fig.~\ref{fig5} is the accuracy obtained using the
state-interaction method, with various basis sizes, for those
cases in which the results in Ref.~\cite{bowmancpl} are
of sufficient precision for comparison. For
the (2,0,0) state the $\{z;x,y\}$ approximation is comparable
in accuracy to a state-interaction calculation with a
39-state basis set consisting of all states
with $n_x+n_y+n_z<5$. For the (1,0,0) state it is less
accurate than a 39-state calculation but more accurate
than a 20-state calculation ($n_x+n_y+n_z<4$).
For the ground state it is less accurate than a 20-state
calculation but more accurate than a 10-state calculation
($n_x+n_y+n_z<3$). For (0,1,0) the $\{x;y,z\}$ approximation
is between the 20- and 39-state calculations in accuracy.
\section{Computational cost}
\label{sec:cost}
One advantage of SCF perturbation theory over exact perturbation
theory is a decrease in computer memory
requirements. Consider a system with $N$ degrees of freedom. For
exact perturbation theory the wavefunction component $\Psi_k$ at
order $k$ in the expansion can be described \cite{dunn} by a
tensor ${\bf a}_k$ of rank $N$ that contains the
coefficients of a linear superposition of harmonic
eigenfunctions. In each dimension of the tensor the number of
nonzero elements increases linearly with $k$. Therefore the SCF
approximation reduces the scaling with $k$ from $k^N$ to
$\sum_\alpha k^{N_\alpha}$ where $N_\alpha$ is the number of
coupled coordinates in coordinate group $\alpha$.
Another advantage of SCF approximations is a decrease in
computational time; that is, the number of arithmetic operations
needed to compute the energy expansion to given asymptotic
order. Consider the computational cost of exact
(unseparated) perturbation theory. The recursion relation
\cite{anharm,dunn} that expresses the tensor ${\bf a}_j$ in terms
of the ${\bf a}_{i