In this communication, we propose a method for obtaining isolated excited states within the full configuration interaction quantum Monte Carlo framework. This method allows for stable sampling with respect to collapse to lower energy states and requires no uncontrolled approximations. In contrast with most previous methods to extract excited state information from quantum Monte Carlo methods, this results from a modification to the underlying propagator, and does not require explicit orthogonalization, analytic continuation, transient estimators, or restriction of the Hilbert space via a trial wavefunction. Furthermore, we show that the propagator can directly yield frequency-domain correlation functions and spectral functions such as the density of states which are difficult to obtain within a traditional quantum Monte Carlo framework. We demonstrate this approach with pilot applications to the neon atom and beryllium dimer.

Almost all of the many variants of projector quantum Monte Carlo (QMC) rely on the properties of the operator e−βH, where due to its similarity to the time-evolution operator, the variable β is denoted “imaginary time.” Generally, this imaginary time is discretized, and the operator iteratively applied as a short-time propagator, in order to simulate its action in the large β limit.1 Expressing an initial wavefunction in the eigenbasis of the Hamiltonian of interest, the application of this e−βH propagator results in a projection onto the

$i{\textrm {th}}$
i th eigenstate proportional to
$e^{-\beta E_i}$
eβEi
, where Ei is the energy of this eigenstate. It is clear to see that in this large β limit, the projection onto the lowest energy eigenvector dominates the wavefunction, whereby ground state properties can be extracted, assuming some overlap with the initial wavefunction. While this formalism is clearly powerful, by construction it exponentially quickly projects out excited states of the system which may be of interest, and are of critical importance in the simulation of finite-temperature properties, reaction dynamics, photochemistry and many other areas.

To date, isolating excited states of systems via projector QMC methods has only been practical with a restriction on the projection to sample a space which is approximately orthogonal to those of the lower energy states, via nodal constraint,2 or orthogonalization against them in a subspace projection method.3,4 More indirectly, statistical methods have been used on short periods of imaginary-time in order to isolate individual decay rates in the spectrum by analytic continuation to a real-time dynamic.5,6 However, these approaches are not entirely satisfactory; accurate nodal surfaces for excited states can be difficult to obtain, resulting in a larger fixed node error and potentially transient estimators, while the subspace projection method has limited applicability.7 In addition, the Bayesian techniques which rely on maximizing entropy to approximate a notoriously unstable inverse Laplace transform, have difficulty achieving quantitative accuracy within noisy data sets.8,9 Despite this, there are examples of accurate excited states within the nodal constraint,10 while maximum entropy techniques are particularly prevalent in solid state calculations to obtain the density of states, often in the case of continuous-time QMC as applied to quantum impurity models and dynamical mean-field theory.11,12

Here, we take a different approach to the calculation of excited states, within the context of full configuration interaction quantum Monte Carlo (FCIQMC).13–15 This recently introduced method applies the imaginary-time evolution propagator to a stochastic “walker” representation of the wavefunction expressed in the full space of Slater determinants constructed from a single-particle basis of size M. Although this reintroduces a basis set error compared to those methods operating in real space, it confers various advantages which mitigate this. The discrete basis allows for an efficient walker annihilation algorithm, which can exactly overcome the fermion sign problem in the sampling, provided enough walkers are used.16 The “initiator” approximation was formulated to maintain a high annihilation rate, and control the growth of noise in a systematically improvable fashion.14,15 This has allowed far larger systems to be treated at an accuracy comparable to that of FCI (exact diagonalization), within small and systematically improvable error bars. In addition, a semi-stochastic adaptation of the algorithm,17 as well the introduction of a partial nodal constraint18 and ideas from quantum chemistry19 hold promise of increased accuracy and efficiency of the method.

In order to project out a targeted excited state rather than the ground state, we propose the use of a projection operator of the form

\begin{equation}P(H)=e^{-\beta ^2(H-S)^2} .\end{equation}
P(H)=eβ2(HS)2.
(1)

For sufficiently large β, this Gaussian propagator will result in the dominant eigenstate being the one closest in energy to the chosen value of the diagonal offset S, termed the shift. In the eigenbasis of H, {|Ψi⟩, Ei} with Ψ0 representing the ground state, and starting from an initial wavefunction

${ |\psi _{\textrm {T}} \rangle }$
|ψT with S = Ek, it can be seen that the long time propagation results in

\begin{equation}|\Psi _k\rangle \propto \lim _{\beta \rightarrow \infty } \sum _i |\Psi _i\rangle e^{-\beta ^2(E_i-E_k)^2} \langle \Psi _i { |\psi _{\textrm {T}} \rangle }\propto \sum _i \delta _{E_i,E_k} |\Psi _i \rangle .\end{equation}
|Ψklimβi|Ψieβ2(EiEk)2Ψi|ψTiδEi,Ek|Ψi.
(2)

We note here that a projector of this form was proposed back in 1983 within continuum QMC approaches,20,21 although due to sign problems, and significantly larger time step errors resulting from the fact that H2 is more singular than H, only one-electron systems were reported, and no modern implementation exists in the literature. This issue of the time step highlights another advantage of working in a finite basis, in that the spectrum is bounded both from below and above. This allows for linearization of the short-time propagator,

\begin{eqnarray}|\Psi \rangle &\propto & \lim _{P\rightarrow \infty } [A e^{-\tau^2 (H-S)^2}]^P {|\psi_{\textrm {T}} \rangle},\end{eqnarray}
|ΨlimP[Aeτ2(HS)2]P|ψT,
(3)
\begin{eqnarray}|\Psi \rangle &\propto & \lim _{P\rightarrow \infty } [A(1-(\tau (H-S))^2)]^P { |\psi _{\textrm {T}} \rangle },\end{eqnarray}
|ΨlimP[A(1(τ(HS))2)]P|ψT,
(4)

without becoming unbound and dominated by very high energy states oscillating in time, and without incurring time step errors in the final wavefunction so long as the time step is less than an upper bound given by

$\tau \le \frac{\sqrt{2}}{E_{\mathrm{max}}-E_{\mathrm{min}}}$
τ2E max E min ⁠. Repeated application of the short-time propagator therefore results in a “power-method” for states on the interior of the spectrum, rather than at the extrema, with similarities to filter diagonalization.22,23

Propagation with Eq. (1) leads to a theoretical decay of state j from i as

$e^{-\beta ^2((E_i-S)^2-(E_j-S)^2)}$
eβ2((EiS)2(EjS)2)⁠. In contrast with the ground state propagator, this rate of decay depends on S, with it being advantageous to choose S to be as close to the energy of the state of interest as possible. However, even if S is chosen exactly, the projection of the non-dominant states is slower compared to the ground state propagation. In addition, unless S is chosen exactly, the long-time propagation of the dynamic will result in a continued projection onto a decaying function of all states, including the dominant one. For this reason, the factor of A is introduced into the short-time propagator, such that at convergence, its value can be varied in order to maintain a constant L1 normalization of the dominant wavefunction and a stable number of walkers. This is analogous to the variation of S in the ground state projection.13 

The differential formulation of the exact dynamic governed by this propagator for a given component of the wavefunction, Ci, can be written as

\begin{equation}\frac{d C_i}{d \tau ^2} = (A-1)C_i - \epsilon A \sum _{j,k} (H_{ij}-\delta _{ij}S)(H_{jk}-\delta _{jk}S) C_k ,\end{equation}
dCidτ2=(A1)CiεAj,k(HijδijS)(HjkδjkS)Ck,
(5)

where ε → τ2 as A → 1, and the application of H2 has been decomposed by a resolution of the identity over the connecting space of determinants j. This formulation is now amenable to stochastic integration with a discrete walker representation of the determinantal wavefunction coefficients C. As with the ground state projection, there is no unique stochastic algorithm for this dynamic, but the one which we found to be efficient involves a double spawning cycle, which requires little additional overhead compared to the ground state algorithm, and no additional memory requirements. Each iteration, the determinants represented by k are run through, and a spawning attempted to determinant j, in the same fashion as the ground state propagation, but in negative time. This results in a spawning probability to a connected determinant j with a stochastically realised signed amplitude of

$\frac{\tau H_{jk}}{P(k|j)}$
τHjkP(k|j)⁠, where P(k|j) represents the normalised probability to randomly select symmetry-connected determinant j from k. Successfully spawned walkers are subsequently propagated again in the same iteration with a now forwards-time signed amplitude of
$-\frac{\tau H_{ji}}{P(j|i)}$
τHjiP(j|i)
. Care must be taken that for determinant k, the diagonal “death” processes from the first application of H are now interpreted as spawning events, which are also subsequently propagated via Hij. Each iteration, the factor of A is applied initially as a separate enhancement of the local population of each determinant, with the absolute population on each determinant growing with probability ACk.

In Fig. 1, we present an illustrative example of the algorithm for the helium dimer in a cc-pVDZ basis, small enough such that the full spectrum of eigenstates can be calculated and the convergence of the method analysed. The value of S was fixed at

$\sim\! 40\, {\textrm {mE}_{\textrm {h}}}$
40 mE h higher than the second excited state, but such that it remained the dominant state in the dynamic, and was subsequently found to be correctly projected out over time. This is despite working in a canonical representation, and starting with a single walker on the Hartree–Fock determinant which had an initial overlap with the ground state of close to one. In order to grow the walkers, A was initially fixed at a value of 1.004, and was varied when a target number of walkers was reached, in common with the procedure for the ground state propagation. The convergence of the projected energy, as defined in Ref. 13, is shown in Fig. 2 for the same simulation, and reflects the decay from the sampled wavefunction of the first excited state.

FIG. 1.

Convergence of the propagation to the second excited state of He2 at 2.5 Å separation in a cc-pVDZ basis. S was fixed at −3.65

${\textrm {E}_{\textrm {h}}}$
Eh⁠, and A at 1.004 until 10 000 walkers were present, denoted by the dotted line, where A was varied in order to keep this number constant. After variation, the average value of A − 1 was 1.4(6) × 10−6.

FIG. 1.

Convergence of the propagation to the second excited state of He2 at 2.5 Å separation in a cc-pVDZ basis. S was fixed at −3.65

${\textrm {E}_{\textrm {h}}}$
Eh⁠, and A at 1.004 until 10 000 walkers were present, denoted by the dotted line, where A was varied in order to keep this number constant. After variation, the average value of A − 1 was 1.4(6) × 10−6.

Close modal
FIG. 2.

Convergence of the projected energy estimate to the exact eigenvalue. The reference determinant for the projection was dynamically adjusted to project onto the largest weighted determinant in the sample.

FIG. 2.

Convergence of the projected energy estimate to the exact eigenvalue. The reference determinant for the projection was dynamically adjusted to project onto the largest weighted determinant in the sample.

Close modal

In order to reliably extend to larger molecular systems, it is worth considering how to transfer the salient elements of the initiator approximation into this new propagation. The basis of this approximation is to attempt to propagate walkers corresponding to wavefunction signal exactly, while walkers judged to be potentially noise are propagated with a truncated Hamiltonian which acts only over the instantaneously occupied subspace.14,15 This is systematically improvable as the instantaneously occupied subspace grows, or the criterion for walkers corresponding to signal becomes more inclusive. The separation between walkers corresponding to signal and potential noise is not unique. However, it seems sensible to retain the tested feature from ground state propagation that newly spawned walkers on previously unoccupied determinants (i) at the end of an iteration, must have come from an initial determinant (k) which is deemed to have a well-established sign, and therefore a population of walkers above

${n_{\textrm {add}}}$
n add ⁠. Consequently, all walkers from the application of the first Hamiltonian operator are kept, while the information as to whether Ci is above the initiator criterion is passed through to the annihilation stage of the final set of spawned walkers. No walkers are therefore aborted over the resolution of the identity between H2 (determinants j in Eq. (5)).

To test this on a larger system, we study an excited state deep in the spectrum of the 10-electron neon atom in a cc-pVDZ basis, with an energy of approximately

$2.5{\textrm {E}_{\textrm {h}}}$
2.5Eh above the ground state. Setting S to equal the CISDTQ energy for the corresponding state, and while remaining in a canonical Hartree–Fock basis and starting from a random distribution of walkers throughout the whole space, we achieved a converged energy of −126.2118(4)
${\textrm {E}_{\textrm {h}}}$
Eh
, compared to the FCI value of −126.21177
${\textrm {E}_{\textrm {h}}}$
Eh
. This value is
$4.76{\textrm {mE}_{\textrm {h}}}$
4.76 mE h
lower than the initial guess provided by CISDTQ. It would be highly advantageous to develop a robust algorithm for varying S dynamically during the run, as is done for the ground state algorithm. This could be used to maximise the rate of growth of walkers or alternatively minimize A, both of which should adjust S to more closely match the eigenvalue of the state, remove reliance on the initial guess and increase the convergence rate. However, since this requires finding a minimum in a quadratically varying and noisy dataset, no robust algorithm has been identified so far.

Dynamic correlation and response functions due to some perturbation, either in the frequency or time domain, are of critical importance in electronic structure theory,24 and are directly measured in experiments to probe the electronic properties of materials through techniques such as neutron scattering or photoelectron spectroscopy. Many methods, including in general QMC approaches, have significant difficulty in calculating these quantities,25 often having to rely on unstable analytic continuation from imaginary time correlation functions,5,6,8,9 while other methods can bias towards high or low energy regimes.12 Although other correlation functions are possible, here we look at the example of an advanced Green's function, a central concept in electronic structure where the “perturbation” at time t = 0 is the creation of a hole in orbital j. For negative time periods, t, these can be written in the time and frequency domain, respectively, as

\begin{eqnarray}G^-(i,j,t) &=& i\langle \Psi _0|a^{\dagger }_i e^{-i(H-E_0-i\delta )t} a_j| \Psi _0 \rangle ,\end{eqnarray}
G(i,j,t)=iΨ0|aiei(HE0iδ)taj|Ψ0,
(6)
\begin{eqnarray}G^-(i,j,\omega ) &=& \langle \Psi _0|a^{\dagger }_i \frac{1}{\omega -(H-E_0)+i\delta } a_j| \Psi _0 \rangle .\end{eqnarray}
G(i,j,ω)=Ψ0|ai1ω(HE0)+iδaj|Ψ0.
(7)

Unlike the inverse Laplace transform required for the analytic continuation of imaginary time correlation functions, the transform between these two domains is a well-conditioned and numerically stable Fourier transform in the presence of noisy data. Spectral density functions, such as the density of states for extended systems, are then defined in the Lehmann representation26 as

\begin{eqnarray}\hspace*{-9.1pc} A^-(i,j,\omega ) &=& -\frac{1}{\pi } \Im [G^-(i,j,\omega )]\end{eqnarray}
A(i,j,ω)=1π[G(i,j,ω)]
(8)
\begin{eqnarray}\hspace*{15pt} &=& \frac{1}{\pi }\sum _{n}\frac{\big\langle \Psi _0^{N}|a^{\dagger }_i|\Psi _n^{N-1}\big\rangle \delta \big\langle \Psi _n^{N-1} | a_j | \Psi _0^{N} \big\rangle }{\big(\omega -E_n^{N-1}+E_0^N\big)^2+\delta ^2} ,\end{eqnarray}
=1πnΨ0N|ai|ΨnN1δΨnN1|aj|Ψ0NωEnN1+E0N2+δ2,
(9)

which in the small δ limit tends to

\begin{align}A^-(i,j,\omega ) =& \sum _{n}\big\langle \Psi _0^{N}|a^{\dagger }_i|\Psi _n^{N-1}\big\rangle \big\langle \Psi _n^{N-1} | a_j | \Psi _0^{N} \big\rangle \nonumber \\&\times\,\delta \big(\omega -(E_n^{N-1}-E_0^N\big)),\end{align}
A(i,j,ω)=nΨ0N|ai|ΨnN1ΨnN1|aj|Ψ0N×δω(EnN1E0N),
(10)

where δ in the above equation represents the dirac-delta function.

Assuming A = 1, application of the propagator in Eq. (1) for a time

$\beta ^2=\frac{1}{2 \delta ^2}$
β2=12δ2 will result in the wavefunction

\begin{equation}C(\beta ^2)=e^{-\frac{1}{2 \delta ^2}(H-S)^2}{ |\psi _{\textrm {T}} \rangle },\end{equation}
C(β2)=e12δ2(HS)2|ψT,
(11)

which when applied to an initial wavefunction

${ |\psi _{\textrm {T}} \rangle }=a_j|\Psi _0 \rangle$
|ψT=aj|Ψ0 obtained from the ground-state dynamic, and then projected onto
$\frac{\beta }{\sqrt{\pi }}\langle \Psi _0 |a^{\dagger }_i$
βπΨ0|ai
will result in the distribution

\begin{align}f(i,j,\omega ) =& \frac{1}{\sqrt{2 \pi } \delta } \sum _{n}\big\langle \Psi _0^{N}|a^{\dagger }_i|\Psi _n^{N-1}\big\rangle \big\langle \Psi _n^{N-1} | a_j | \Psi _0^{N} \big\rangle \nonumber \\& \times\,e^{-\frac{1}{2 \delta ^2}(\omega -E_n^{N-1}+E_0^N)^2,}\end{align}
f(i,j,ω)=12πδnΨ0N|ai|ΨnN1ΨnN1|aj|Ψ0N×e12δ2(ωEnN1+E0N)2,
(12)

for

$S=\omega +E_0^N$
S=ω+E0N⁠. This will tend to the spectral function given in Eq. (10) in the large β limit. The real parts of the Green's function can then be obtained if needed from the Kramers-Kronig relation.24 Results from a pilot investigation of the beryllium dimer in a cc-pVDZ basis, where the exact Green's function can be obtained from complete diagonalization, are shown in Fig. 3.

FIG. 3.

High energy window of the spectral function A(1, 1, ω) as given in Eq. (9) with

$\delta =0.0141{\textrm {E}_{\textrm {h}}}$
δ=0.0141Eh⁠, and stochastic evaluation via FCIQMC for an equivalent time β = 50 a.u. for frozen-core Be2 in a cc-pVDZ basis at 2.5 Å. Vertical lines indicate the difference between the ground state energy and the eigenvalues of the N-1 system symmetry connected in G(1, 1, ω), although some are coupled too weakly to contribute significantly to the spectral function. Approximately 10 independent calculations at each value of ω were averaged to obtain the errorbars.

FIG. 3.

High energy window of the spectral function A(1, 1, ω) as given in Eq. (9) with

$\delta =0.0141{\textrm {E}_{\textrm {h}}}$
δ=0.0141Eh⁠, and stochastic evaluation via FCIQMC for an equivalent time β = 50 a.u. for frozen-core Be2 in a cc-pVDZ basis at 2.5 Å. Vertical lines indicate the difference between the ground state energy and the eigenvalues of the N-1 system symmetry connected in G(1, 1, ω), although some are coupled too weakly to contribute significantly to the spectral function. Approximately 10 independent calculations at each value of ω were averaged to obtain the errorbars.

Close modal

In order to reduce the statistical error, it may be necessary to average over a small number of independent calculations at each frequency, and this can be combined with an elimination of the bias derived from choosing a correlated sample of

$\langle \Psi _0|a^{\dagger }_i$
Ψ0|ai and aj0⟩,3 by taking the Ψ0 samples on each side of Eq. (12) from different snapshots in imaginary time. In addition, by storing multiple wavefunctions of the type
$\langle \Psi _0|a^{\dagger }_i$
Ψ0|ai
at the same time, all M2 single-particle Green's functions can be calculated at a cost of
$\mathcal {O}[M]$
O[M]
FCIQMC calculations per frequency point, without the expectation of any variation in accuracy between high and low energy regimes.

However, despite modest successes, it is clear that obtaining converged results through the use of this operator is substantially more difficult than with the ground state projection. This is mainly due to an additional factor of (τΔE)−1 in the number of iterations required to project out undesired states with energy gap ΔE for comparable accuracy to the ground state propagation. The result is that while in the ground state propagation excited states were filtered relatively quickly with only isolated convergence issues in the case of near degeneracy,15 the number of iterations required for excited state propagation are substantially increased, as well as the dynamic being less well-conditioned with respect to walker fluctuations. This is also exacerbated by a generally more multiconfigurational wavefunction which increases random error in the projected energy estimator.13 A more judicious choice of orbital basis and initial conditions optimized for the state of interest, as well as a multireference projected energy formulation17 would ameliorate many of these issues. In addition, there is the possibility of preconditioning techniques familiar from iterative diagonalization methods27 being transferred into the stochastic dynamic, as well other operators, such as e−β|H|, which should behave more efficiently and allow for extension to larger systems. Research in these directions is currently under way. It is clear that alternative propagators within the FCIQMC dynamic holds promise for obtaining accurate excited states.

We would like to acknowledge financial support from the U.S. Department of Energy, Office of Science, through Award Nos. DE-SC0008624 and DE-FG02-07ER46432. We thank Cyrus Umrigar for helpful comments.

1.
W. M. C.
Foulkes
,
L.
Mitas
,
R. J.
Needs
, and
G.
Rajagopal
,
Rev. Mod. Phys.
73
,
33
(
2001
).
2.
R.
Grimes
,
B.
Hammond
,
P.
Reynolds
, and
W. J.
Lester
,
J. Chem. Phys.
85
,
4749
(
1986
).
3.
D.
Ceperley
and
B.
Bernu
,
J. Chem. Phys.
89
,
6316
(
1988
).
4.
Y.
Ohtsuka
and
S.
Nagase
,
Chem. Phys. Lett.
485
,
367
(
2010
).
5.
G.
Baym
and
N.
Mermin
,
J. Math. Phys.
2
,
232
(
1961
).
6.
J.
Gubernatis
,
M.
Jarrell
,
R.
Silver
, and
D.
Sivia
,
Phys. Rev. B
44
,
6011
(
1991
).
7.
M.
Jones
,
G.
Ortiz
, and
D.
Ceperley
,
Phys. Rev. E
55
,
6202
(
1997
).
8.
D.
Thirumalai
and
B.
Berne
,
J. Chem. Phys.
79
,
5029
(
1983
).
9.
O.
Gunnarsson
,
M. W.
Haverkort
, and
G.
Sangiovanni
,
Phys. Rev. B
81
,
155107
(
2010
).
10.
P.
Zimmerman
,
J.
Toulouse
,
Z.
Zhang
,
C.
Musgrave
, and
C.
Umrigar
,
J. Chem. Phys.
131
,
124103
(
2009
).
11.
Y.-H.
Chen
,
H.-S.
Tao
,
D.-X.
Yao
, and
W.-M.
Liu
,
Phys. Rev. Lett.
108
,
246402
(
2012
).
12.
E.
Gull
,
A.
Millis
,
A.
Lichtenstein
,
A.
Rubtsov
,
M.
Troyer
, and
P.
Werner
,
Rev. Mod. Phys.
83
,
349
(
2011
).
13.
G. H.
Booth
,
A. J. W.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
131
,
054106
(
2009
).
14.
D.
Cleland
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
132
,
041103
(
2010
).
15.
G. H.
Booth
,
D.
Cleland
,
A. J. W.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
135
,
084104
(
2011
).
16.
J. S.
Spencer
,
N. S.
Blunt
, and
W. M. C.
Foulkes
,
J. Chem. Phys.
136
,
054110
(
2012
).
17.
F.
Petruzielo
,
A.
Holmes
,
H.
Changlani
,
M.
Nightingale
, and
C.
Umrigar
, “
Semistochastic projector Monte Carlo method
,”
Phys. Rev. Lett.
(to be published); preprint arXiv:1207.6138.
18.
M.
Kolodrubetz
and
B.
Clark
,
Phys. Rev. B
86
,
075109
(
2012
).
19.
G. H.
Booth
,
D.
Cleland
,
A.
Alavi
, and
D. P.
Tew
,
J. Chem. Phys.
137
,
164112
(
2012
).
20.
J.
Hirsch
and
J.
Schrieffer
,
Phys. Rev. B
28
,
5353
(
1983
).
21.
D.
Thirumalai
,
B. C.
Garrett
, and
B. J.
Berne
,
J. Chem. Phys.
83
,
2972
(
1985
).
22.
G.
Grosso
,
L.
Martinelli
, and
G.
Parravicini
,
Nuovo Cimento D
15
,
269
(
1993
).
23.
M.
Wall
and
D.
Neuhauser
,
J. Chem. Phys.
102
,
8011
(
1995
).
24.
E.
Gagliano
and
C.
Balseiro
,
Phys. Rev. Lett.
59
,
2999
(
1987
).
25.
D.
Thirumalai
and
B.
Berne
,
Comput. Phys. Commun.
63
,
415
(
1991
).
26.
L.
Fetter
and
J.
Walecka
, in
Quantum Theory of Many-Particle Systems
(
Dover
,
2003
).
27.
E.
Davidson
,
J. Comput. Phys.
17
,
87
(
1975
).