Nonspectral Methods for Solving the Schrödinger Equation for

Aug 5, 2011 - This Perspective briefly presents nonspectral (pseudospectral and grid) methods to solve the Schrödinger equation and their recent and ...
0 downloads 0 Views 2MB Size
Download PDF
PERSPECTIVE pubs.acs.org/JPCL

€ dinger Equation for Nonspectral Methods for Solving the Schro Electronic and Vibrational Problems Sergei Manzhos,*,† Tucker Carrington, Jr,*,‡ and Koichi Yamashita*,§ †

Research Center for Advanced Science and Technology (RCAST), The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8904, Japan ‡ Chemistry Department, Queen’s University, Kingston, Ontario, K7L 3N6, Canada § Department of Chemical System Engineering, School of Engineering, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan ABSTRACT: To compute molecular and material properties, one must solve the electronic and/or vibrational Schr€odinger equation. This Perspective briefly presents nonspectral (pseudospectral and grid) methods to solve the Schr€odinger equation and their recent and promising applications to vibrational and electronic problems. Pseudospectral and grid methods facilitate parallelization, which is of critical importance for modern material modeling applications. By obviating the need to choose basis functions for which exact integrals are possible, they also enable one to use smaller basis sets. For the vibrational problem, pseudospectral and grid methods, in principle, make it possible to compute a vibrational spectrum without a potential surface, simply from values of the potential at points. This is possible if optimized basis functions are used. Iterative eigensolvers are very advantageous when simpler basis functions are used because matrixvector products can be evaluated efficiently.

B

y solving the Schr€odinger equation (SE), chemists are able to predict and understand molecular and material properties. From the solution of the electronic SE (ESE), one obtains potentials, information about reaction pathways, harmonic frequencies, and electronic, optical, magnetic, and mechanical properties of molecules and materials. From the solution of the vibrational SE (VSE), one can obtain infrared (IR) spectra, which make it possible to identify species and compute properties of transition states.

€dinger equation By solving the Schro (SE), chemists are able to predict and understand molecular and material properties. A vibrational spectrum is computed by solving the SE ^ nuc þ V^ nuc Þψnuc ðRÞ ¼ Eψnuc ðRÞ ^ nuc ψnuc ðRÞ  ðT H

ð1Þ

^ nuc is a Hamiltonian operator, which consists of a nuclear where H ^ nuc (KEO) and a potential energy kinetic energy operator T surface (PES), and V^ nuc = Vnuc(R), R = (R1, R2, ..., Rd). Here ,d is the number of vibrational coordinates; d = 3I  6 for an isolated molecule with I atoms, and d = 3I for a molecule in an external potential (such as a molecule on a surface). According to the BornOppenheimer approximation, Vnuc(R) is a solution of the r 2011 American Chemical Society

electronic Schr€odinger equation ^e þ V ^ el ψel ðrÞ ¼ ðT ^ en þ U ^ ee Þψel ðrÞ H ¼

N 1 N 2 ∇ þ Ven ðri Þ þ  2 i¼1 i i¼1

∑

∑

N

!

∑ Uee ðri , rj Þ i1 |rj  ri|1 is the Coulomb repulsion between the electrons. We use atomic units unless stated otherwise. The BornOppenheimer approximation enables one to divide the full problem into two pieces. It is a good starting point even when nonadiabatic effects1 are important. In this Perspective, we discuss and compare methods for solving the VSE and the ESE. Methods for solving the VSE and the ESE have evolved along somewhat different tracks for several reasons. (1) To solve the ESE, it is possible to choose good basis (Gaussian-based) functions with which all integrals can be calculated exactly. To solve the VSE, numerical evaluation of integrals is imperative. (2) Whereas when solving the VSE, one usually wishes for hundreds of energies, when solving the ESE, one requires one or only a few electronic energy levels. (3) The number of vibrational degrees of freedom (DOF) is much less Received: April 15, 2011 Accepted: August 5, 2011 Published: August 05, 2011 2193

dx.doi.org/10.1021/jz200513h | J. Phys. Chem. Lett. 2011, 2, 2193–2199

The Journal of Physical Chemistry Letters

PERSPECTIVE

2 Here, F(r) = ∑N i=1 |ϕi(r)| is the density, Ven(r) has the same meaning as in eq 2, and the exchange-correlation potential Vxc[F(r)] includes many-particle effects. As the last two terms depend on ϕi(r) in a nonlinear way, the solution proceeds by selfconsistent iterations. The introduction of PP greatly simplifies the solution of the ESE and makes possible calculations on systems with thousands of atoms and therefore computational material design. In the pseudopotential approximation (or the related projector augmented-wave (PAW)11 scheme), only valence electrons are treated explicitly, and the effect of the core is modeled by replacing the (singular) nuclear potential with a well-behaved function12 (see the abstract graphic). The valence electrons’ wave functions are thus made smooth enough to be expanded in a simple basis. When DFT and PP are used, eqs 3 and 1 with Veff(r) and Vnuc(R) are similar. Consequently, comparing methods developed to solve them is useful. ^ The most common approach for solving Hψ(x) = Eψ(x) is to convert it to a matrix problem by substituting a basis expansion of the wave function

ψðxÞ ¼

M

∑ ci Φi ðxÞ i¼1

ðG þ V RÞc ¼ RcE0

)

upper bounds for energy levels. As the basis size M is increased to minimize the expansion error ε = ψ(x)  ∑M i=1 ciΦi(x) , the eigenvalues E0 converge to E; moreover, an error ε in a wave function corresponds to an error ε2 in the energy.13 When one abandons the idea of using exact matrix elements, it becomes possible to choose from a large variety of bases. In vibrational calculations, it is common to use 1-d bases, each of whose functions is the product of a polynomial and the square root of a weight function.14 In such a basis, overlap and KEO matrix elements are known exactly, but potential matrix elements are not. If a PP is used, plane wave basis functions have important advantages for electronic calculations.57,9,10 Kinetic and overlap matrix elements are exact, and for extended systems, periodic boundary conditions (PBC) are easily imposed. They also avoid basis set superposition errors (BSSEs) and obviate the need to consider Pulay terms. The main advantage of plane wave functions is their simplicity; a single parameter controls the basis size. For both the VSE and the ESE, simple basis functions have the disadvantage that the size of the basis required to obtain converged results is large. Using atom-centered localized basis sets, the ESE can be solved with direct diagonalization for systems with as many as hundreds of atoms (M ≈ 104), but when plane wave functions are used, for bigger systems, it becomes necessary to use an iterative eigensolver.15,16 If products of polynomial basis functions are used, iterative eigensolvers are necessary for any molecule with more than three atoms. For information about iterative eigensolvers, see refs 14 and 1719. For both DFT and the VSE, it is not possible, regardless of the choice of the basis, to use exact matrix elements. There are several alternative methods, collocation, pseudospectral (PS), neural network (NN), finite difference (FD), and so forth.16 All of these have advantages. Given the wide and increasing availability of huge parallel computers, parallelizability is an important property of a good method. The parallelizability of some codes that do not use FD is poor. This can be an important disadvantage even when only dozens of CPUs are used. For example, Figure 1 shows the scaling that can be achieved with a popular electronic structure code after careful selection of parameters. All of the nonspectral methods can be parallelized. In the collocation approach,2124 the wave function is expanded as in eq 4 and substituted into the SE, but rather than multiplying by a basis function on the left, one multiplies by a delta function and then integrates. The matrix eigenvalue problem is )

than the number of electronic DOF. (4) The vibrational KEO is more complicated because coordinates are used that enable one to exploit the conservation of angular momentum. (5) Because the self-consistent field approach works so well, it is almost always used, and therefore, the ESE that one usually solves is actually nonlinear. (6) The potential in the ESE has multiple singularities. (7) Frequently, one wishes to solve the ESE for many molecular shapes. (8) Frequently, the absolute error considered acceptable is larger for the ESE than that for the VSE. For the VSE, an acceptable error is usually less than 1 cm1; for the ESE, it is usually larger. Although there are still important differences, lower-accuracy electronic structure methods and VSE methods are fairly similar. If density functional theory (DFT) and pseudopotentials (PP) are used, this is especially true. DFT has become the method of choice for systems with hundreds or thousands of atoms and, unless very accurate calculations are required, is also often used for systems with fewer atoms. All codes used by material scientists employ DFT and PP.210 In DFT, rather than solving a many-electron SE, eq 2, one solves a single-particle (3-d) KohnSham (KS) equation with an effective potential   ^ KS ϕi ðrÞ ¼ 1∇2 þ Veff ðrÞ ϕi ðrÞ ¼ εi ϕi ðrÞ H 2 Z ð3Þ Veff ðrÞ ¼ Ven ðrÞ þ dr 0 Fðr 0 Þjr  r 0 j1 þ Vxc ½FðrÞ

ð6Þ

ð4Þ

into the SE, multiplying on the left with one of the basis functions used to expand ψ(x) and integrating. This yields the matrix equation Hc ¼ E0 Sc

ð5Þ

^ jæ, the Hamiltonian matrix, and Sij = ÆΦi|Φjæ, with Hij = ÆΦi|H|Φ the overlap matrix. If matrix elements are exact, we call this a spectral method. When it is possible to choose good basis functions in terms of which there are analytic equations for the matrix elements, the spectral approach is a good option. For small molecules, this method of solving the ESE has been used with great success. It has the advantage that eigenvalues of eq 5 are

Figure 1. Scaling of a SIESTA calculation for a pilin protein of 944 atoms with direct diagonalization and a DZP basis, as reported in ref 20. 2194

dx.doi.org/10.1021/jz200513h |J. Phys. Chem. Lett. 2011, 2, 2193–2199

The Journal of Physical Chemistry Letters

PERSPECTIVE

^ Φj(xi), and Vij = V(xi)δij. To build the where Rij = Φj(xi), Gij = T matrices, one needs only to evaluate V(x) at points xj; no integrals are performed. Multiplying eq 6 on the left by RTW, where W is a (any) diagonal weight matrix, one sees that the collocation eigenvalues are equal to the eigenvalues that one would obtain by expanding as in eq 4, substituting into the SE, multiplying by a Φi(x) basis function on the left, integrating, and doing all integrals by quadrature. The error in the jth collocation energy level is therefore proportional to ε2j + εqεj, where εq is a measure of the quadrature error.25 If the collocation points are chosen to make the quadrature error small, the error in the collocation energies is nearly ε2j ; if the points are poorly chosen, the error may be much larger. To use collocation, it is necessary to have as many points as basis functions. Therefore, it is not possible to take full advantage of the small size of a good basis when an equally small but accurate quadrature point set is not available. Because it is hard to devise good quadratures, one is often forced to use a dense grid of points and therefore a large basis. See, however, refs 26 and 27. The most serious problem with collocation is the fact that the associated matrix eigenvalue problem is both generalized and nonsymmetric. This makes its solution costly (less so if R is diagonal). The collocation method has been used to solve the VSE but only for 2- and 3-d problems.21,22,26 Quadrature of some kind is a popular alternative to collocation. The quadrature eigenvalue problem obtained from collocation is ðRT W G þ R T W V RÞc ¼ RT W RcE0

ð7Þ

Frequently, RTWG and RTWR are replaced with exact matrices, and only the potential matrix elements are obtained with quadrature. This is known as the finite basis representation (FBR).14 The application of RTWVR to a vector can be thought of as proceeding in steps: transform the vector to the grid, apply V, and transform28 back to the basis. This basic pseudospectral idea is used in many forms. For vibrational problems, it is common to evaluate potential integrals with product Gauss quadrature. This makes it possible to transform from the basis to the grid coordinate-by-coordinate and to use sequential summation.2932 A crucial advantage of the FBR is the fact that matrix vector products are done in sequential steps (rather than computing potential matrix elements and then multiplying the potential matrix by the vector). Because the product quadrature grid is large, the potential is needed at many points. It is therefore crucial to reduce the number of quadrature points.26,27 To solve the ESE, it is the effective potential that must be computed on the grid. Matrices in the discrete variable representation14 (DVR) are unitarily equivalent to FBR matrices, and the DVR can also be thought of as a pseudospectral method. In electronic structure theory, Friesner3439 showed that a PS approach similar to eq 7 is about as accurate as the traditional spectral method. Rather than multiplying the collocation equation by RTW, it is multiplied by R1. When used with the HartreeFock method, the most important advantage of the PS technique is due to the fact that two-electron integrals are not computed explicitly (cf. that in the FBR matrix, elements of the potential are not computed). Although this is not relevant for DFT, the PS approach is still computationally efficient for essentially the same reasons as those for the vibrational problem. It is possible to both vectorize and parallelize PS calculations.40 Similar ideas have been applied to methods that account for

Figure 2. Parallel speedup of different parts of a pseudospectral algorithm of Chasman et al. for an alanine pentapeptide calculation reported in ref 33. The transformation between the spectral and real space (Least squares) is not parallelized well. Reprinted with permission from ref 33. Copyright (1998) by John Wiley & Sons, Inc.

correlation,36,4146 demonstrating that properly parallelized PS methods work well.33,40,41,47 Applied to DFT and time-dependent DFT (TDDFT), the PS ideas also resulted in a significant speedup.48,49 Unless the grid and the basis functions are carefully chosen, it is necessary to invert R, which can be costly. Avoiding aliasing errors39 can also be tricky. Both the inversion of R and the dealiasing procedure do not parallelize easily (Figure 2). There are DFT methods that, unlike collocation and PS methods, do not expand the wave function in a basis. Becke’s NUMOL is one example.50 It builds a wave function from solutions to local SEs, each with a spherically averaged potential, that are computed using a spherical harmonic basis and a sevenpoint finite difference approximation for the radial coordinate. One local SE is associated with each atom, and each has its own and radial and angular grid. The method was successfully tested on molecules with up to seven atoms.51 However, for each atom, one must self-consistently solve a local SE, and each has a different spherical harmonic expansion. In addition, one must interpolate local solutions onto each other’s grids. This method is cumbersome for large systems. Simple finite difference for the full problem is also an option.52 For the VSE, FD approaches require many more points than do PS methods. Without PP, FD approaches for the ESE are known to be poor, but with PP, even an equally spaced FD scheme works rather well and has the advantage that it is relatively easy to parallelize.11,53 The required grid spacing and the plane wave cutoff are linked. Adaptive grids can be achieved via a coordinate transformation.54 Simple central finite difference schemes for the KEO have the form ∇2 ψðxi Þ≈

K

∑j ci ψðxjÞ

ð8Þ

where xj includes xi and its K  1 nearest neighbors. This scheme is used in OCTOPUS,55 for example. Higher accuracy can be achieved by using more local information and the Mehrstellen operator53 to discretize the KS equation 1 HMehr ½j ¼ AMehr ½j þ BMehr ½V j ¼ εBMehr ½j 2 2195

ð9Þ

dx.doi.org/10.1021/jz200513h |J. Phys. Chem. Lett. 2011, 2, 2193–2199

The Journal of Physical Chemistry Letters

PERSPECTIVE

where, typically, σ(x) = (1 + ex)1. σ depends on a linear combination of the coordinates but is a univariate function. The parameters w, b, and c are adjusted to minimize the residual M

R ¼

∑ jðH^  EÞjðxi Þj2 i¼1 M

∑ jjðxi Þj

ð11Þ

2

i¼1

Figure 3. Parallel scaling for a ground-state total energy calculation of a 762-atom system consisting of a 102-atom Au cluster surrounded by 44 p-MBA. Reprinted with permission from ref 11. Copyright (2010) by IOP Publishing Ltd. DOI 10.1088/0953-8984/22/25/253202.

where both AMehr and BMehr have the form of eq 8. The Mehstellen stencil is used in GPAW.11 Application of operators in eq 8 or 9 results in sparser matrices and, more importantly, matrices with smaller bandwidths than those in a spectral or PS approaches.56 This permits an effective parallelization of the matrixvector product required to use iterative eigensolvers. As a result, near-linear scaling can be achieved on up to 104 CPUs for systems with more than 700 atoms, as is shown in Figure 3 (cf. Figure 1). Using a grid spacing that is not small enough introduces errors in energies and thereby forces. Even the energy of a single atom will vary as it is moved across a flat potential. This is mitigated by Fourier filtering the potentials20,53 with a plane wave cutoff corresponding to the grid and by averaging over shifted positions of the grid points.20

A plane wave basis necessitates the use of a multi-d fast Fourier transform (FFT), which can be effectively implemented on a single processor but is difficult to parallelize. The difficulty arises from the fact that vector components over which one sums may be located on different processors. The matrixvector product with a banded matrix is parallelized by storing different parts of the matrix and the vector on different processors, and the more banded the matrix, the better the parallelization. Although, as yet relatively untested, neural networks (NN)57 can also be used to compute energy levels. Like other approaches discussed in this perspective, the NN method requires only values of the potential at points. The potential advantage of the NN method is the fact that only a small number of points are required. The wave function is represented as a sum of basis functions σ called “neurons” N

∑ cq σðyq ðbq , wq , xÞÞ þ bÞ q¼1

jk ðxÞ ¼

N

N

∑ cki gi ðb, w, xÞ ¼ i∑¼ 1 cki i¼1

ð10Þ

d Y 2 bij 2 pffiffiffi ebij ðwij  xj Þ π j¼1

ð12Þ As in collocation, we require that the SE be satisfied at a set of points xi, n = 1, ..., M, randomly chosen in the region where the wave function is expected to be non-negligible. Unlike collocation, we choose more points than basis functions. This is necessary because the (optimized) basis set is small. The matrix equation that we obtain is rectangular. ðM  Ek SÞck ¼ 0

Pseudospectral and grid methods facilitate parallelization, which is of critical importance for modern material modeling applications. Near-linear scaling can be achieved on up to 104 CPUs.

jðxÞ ¼ σ 0 ð

The KEO can be applied to the wave function in eq 10 analytically. E can be computed from the wave function58 or adjusted as a fitting parameter.59,60 The strength of this approach is the adaptability of the basis, which significantly reduces its size. The required nonlinear optimization of hundreds of NN parameters for every wave function may, however, be very costly. NNs have so far been used to find several low-lying vibrational levels on lowdimensional model potentials.5861 Recently, an application to the KS equation for light atoms was reported with modest accuracy.62 We developed a grid-based method to optimize R.61,63,64 We used radial basis function neurons, that is, a Gaussian-like basis

ð13Þ

where M and S are rectangular matrices (of dimension M  N) depending on values of the potential and of the basis functions at different points in the configuration space. See refs 61, 63, and 64 for more detail. For given values of b and w, expansion coefficients are elements of the “pseudoeigenvector” of the rectangular eigenvalue problem corresponding to the smallest eigenvalue. The parameters b and w are then optimized. A pseudoeigenvector is an eigenvector of (M  EkS)H(M  EkS). Energies are not approximated by eigenvalues of the rectangular problem, but Ek for each level is a fitting parameter that is adjusted simultaneously with w and b. A key advantage of this approximation is that the size of the (square) matrix eigenvalue problem is the size of the basis, which can be quite small. This NN algorithm was successfully applied to compute vibrational levels of H2O.64 Using a normal coordinate KEO, Manzhos et al. computed vibrational frequencies of H2O on Pt(111)63 (Figure 4). Spectra were computed directly from 1500 values of Vn by ab initio methods (no PES function was required) with only a few dozen basis functions. The method has not yet been applied to electronic problems. An earlier review of collocation, NN, and related methods is presented in ref 65. Nakatsuji has also proposed a grid-based method with rectangular matrices to solve the free complement local Schr€odinger equation (FC LSE).66 A rectangular matrix equation similar to eq 13 was obtained but then was made square to solve a conventional generalized eigenvalue problem for energies and wave function coefficients. The electronic SE for several atoms and diatomic molecules was solved with high accuracy.66,67 2196

dx.doi.org/10.1021/jz200513h |J. Phys. Chem. Lett. 2011, 2, 2193–2199

The Journal of Physical Chemistry Letters

PERSPECTIVE

To use the methods that we have reviewed, it is necessary to store in memory a vector with as many elements as the number of collocation (or quadrature) points. For example, in vibrational spectroscopy, to calculate (numerically exact) energy levels of a sixatom molecule like CH3CN with a product grid, one would need about 3  1013 potential points. For difficult problems, this clearly limits what is computationally possible. One way to deal with this problem is to develop better quadrature schemes.26,27 For CH3CN, it is possible to reduce the number of points by 6 orders of magnitude. Another is to exploit Boys’s observation that if the basis set is excellent, quadrature (or collocation) error is unimportant.25 One option is to combine an adaptable NN basis and a rectangular collocation matrix equation.61,63,65 The ultimate goal is to reduce the number of points to the extent that they can all be calculated by ab initio methods, obviating the need for a PES. If only modest accuracy is required, this is possible.63 This method could also be applied to solve the KS equation, and it could be used with the FC LSE approach.66,67 Parallelizability is a crucial issue. FD schemes that are useful for DFT calculations with PP facilitate parallelization.

’ AUTHOR INFORMATION Figure 4. Anharmonic frequencies of H2O on Pt(111) were computed in ref 63 directly from DFT data with an NN+collocation method. Atoms: blue, Pt; red, O; gray, H. The blue line with circles is a DFT scan of the PES along the symmetric stretch mode. The corresponding harmonic curve is shown as a dotted black line. The green line is the anharmonic wave function for v = 1.

To use DFT, one often solves a Poisson equation to obtain an effective potential. When discretizing the Poisson equation, it may also be advantageous to use a grid rather than a plane wave basis and the FFT. This allows the choice of Dirichlet rather than PBC in any direction. For example, in the newest version of SIESTA, the FFT approach is replaced with a multigrid solver20 that uses a sequence of grids with multiple resolutions. In that way, convergence is improved at little extra cost. This is because on a given grid, frequency components of the error corresponding to the grid resolution are reduced most effectively.53 The multigrid approach was implemented already in the PS method of Friesner35 and is used in the latest real space codes.11,54 In conclusion, similar nonspectral methods used to solve the VSE and the ESE now enable one to solve difficult and important problems. In electronic structure theory, pseudospectral and finite difference methods have improved speed and parallelizability. Pseudospectral methods are used routinely to solve the VSE. The advantages of sequential summation14,29,30,68 and various parallelization strategies27,69,70 are widely appreciated. The methods in both fields are similar, and further development will profit from cross fertilization.

For the vibrational problem, pseudospectral and grid methods, in principle, make it possible to compute a vibrational spectrum without a potential surface, simply from values of the potential at points.

Corresponding Authors

*E-mail: [email protected] (S.M.); Tucker.Carrington@ queensu.ca (T.C.); [email protected] (K.Y.).

’ BIOGRAPHIES Sergei Manzhos received a M.Sc. from Kharkiv National University in 1999 and a Ph.D. from Queen’s University, Canada, in 2004. In 20052007, he was a postdoctoral fellow at the University of Montreal. In 20082009, he was Assistant Professor at The University of Tokyo, Japan, and since 2010, he has been Assistant Professor at RCAST, The University of Tokyo. Research interests include theoretical spectroscopy and theory of moleculesurface reactions and of dye-sensitized solar cells. Tucker Carrington, Jr. is Canada Research Chair for computational quantum dynamics at Queen’s University. Working with Prof. W. H. Miller, he obtained his Ph.D. in 1985 from UC Berkeley. After postdoctoral work, he joined the Chemistry Department at the University of Montreal in 1988. Most of his research is related to the development of new methods for studying the dynamics of molecules and reactions. For more detail, see http://www.chem.queensu.ca/people/faculty/ Carrington/index.asp. Koichi Yamashita has been Professor at the University of Tokyo since 1997. His undergraduate and graduate degrees were from Kyoto University, supervised by Prof. Kenichi Fukui. He was a postdoctoral fellow with Prof. W. H. Miller at the University of California, Berkeley, from 1982 to 1984. Research interests include ab initio methods, quantum reaction dynamics, and modeling of catalysis. For more detail, see http://www.tcl.t. u-tokyo.ac.jp/english_index.html. ’ REFERENCES (1) Fischer, S. A.; Habenicht, B. F.; Madrid, A. B.; Duncan, W. R.; Prezhdo, O. V. Regarding the Validity of the Time-Dependent Kohn Sham Approach for ElectronNuclear Dynamics via Trajectory Surface Hopping. J. Chem. Phys. 2011, 134, 024102. (2) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09; Gaussian, Inc.: Wallingford CT, 2009. 2197

dx.doi.org/10.1021/jz200513h |J. Phys. Chem. Lett. 2011, 2, 2193–2199

The Journal of Physical Chemistry Letters (3) Soler, J. M.; Artacho, E.; Dale, J. D.; Garcia, A.; Junquera, J.; Ordejon, P.; Sanchez-Portal, D. The SIESTA Method for Ab Initio OrderN Materials Simulation. J. Phys.: Condens. Matter 2002, 14, 2745–2779. (4) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. A.; Zicovich-Wilson, C. M. CRYSTAL: A Computational Tool for the Ab Initio Study of the Electronic Properties of Crystals. Z. Kristallogr. 2005, 220, 571–573. (5) Gonze, X.; Amadon, B.; Anglade, P. M.; Beuken, J. M.; Bottin, F.; Boulanger, P.; Bruneval, F.; Caliste, D.; Caracas, R.; Cote, M.; et al. ABINIT: First-Principles Approach to Material and Nanosystem Properties. Comput. Phys. Commun. 2009, 180, 2582–2615. (6) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21. (7) Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. J.; Refson, K.; Payne, M. C. First Principles Methods Using CASTEP. Z. Kristallogr. 2005, 220, 567–570. (8) SCM, Theoretical Chemistry, ADF2009.01; Vrije Universiteit: Amsterdam, The Netherlands, 2009. (9) Kresse, G.; Furthmuller, J. Efficient Iterative Schemes for ab initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. (10) CPMD; http://www.cpmd.org/, Copyright IBM Corp 1990 2008, Copyright MPI f€ur Festk€orperforschung Stuttgart 19972001; (2008). (11) Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dulak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A.; et al. Electronic Structure Calculations with GPAW: A Real-Space Implementation of the Projector Augmented-Wave Method. J. Phys.: Condens. Matter 2010, 22, 253202. (12) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for Plane-Wave Calculations. Phys. Rev. B 1991, 43, 1993–2006. (13) Stewart, G. W. Matrix Algorithms Vol. II: Eigensystems; SIAM: Philadelphia, PA, 2001. (14) Light, J. C.; Carrington, T., Jr. Discrete Variable Represenrtations and their Utilization. Adv. Chem. Phys. 2000, 114, 263–310. (15) Hutter, J.; Luthi, H. P.; Parrinello, M. Electronic Structure Optimization in Plane-Wave-Based Density Functional Calculations by Direct Inversion in the Iterative Subspace. J. Comput. Mater. Sci. 1994, 2, 244. (16) Torsti, T.; Eirola, T.; Enkovaara, J.; Hakala, T.; Havu, P.; Havu, V.; Hoynalanmaa, T.; Ignatius, J.; Lyly, M.; Makkonen, I.; et al. Three Real-Space Discretization Techniques in Electronic Structure Calculations. Phys. Status Solidi B 2006, 243, 1016–1053. (17) Goedecker, S. Linear Scaling Electronic Structure Methods. Rev. Mod. Phys. 1999, 71, 1085–1123. (18) Carrington, T., Jr. Methods for Calculating Vibrational Energy Levels. In Encyclopedia of Computational Chemistry, 2nd ed.; von Rague Schleyer, P., Ed.; John Wiley & Sons: New York, 2003. (19) Davidson, E. R. Super-Matrix Methods. Comput. Phys. Commun. 1989, 53, 49–60. (20) Artacho, E.; Anglada, E.; Dieguez, O.; Gale, J. D.; Garcia, A.; Junquera, J.; Martin, R. M.; Ordejon, P.; Pruneda, J. M.; Sanchez-Portal, D.; Soler, J. M. The Siesta Method; Developments and Applicability. J. Phys.: Condens. Matter 2008, 20, 064208. (21) Peet, A. C.; Yang, W. The Collocation Method for Calculating Vibrational Bound States of Molecular Systems with Application to Ar-HCl. J. Chem. Phys. 1989, 90, 1746–1751. (22) Peet, A. C.; Yang, W. An Adapted Form of the Collocation Method for Calculating Energy Levels of Rotating AtomDiatom Complexes. J. Chem. Phys. 1989, 91, 6598–6603. (23) Yang, W.; Peet, A. C.; Miller, W. H. A Collocation Approach for Quantum Scattering Based on the S-Matrix Version of the Kohn Variational Principle. J. Chem. Phys. 1989, 91, 7537–7542. (24) Yang, W.; Peet, A. C. The Collocation Method for Bound Solutions of the Schroedinger Equation. Chem. Phys. Lett. 1988, 153, 98–104.

PERSPECTIVE

(25) Boys, S. F. Some Bilinear Convergence Characteristics of the Solutions of Dissymmetric Secular Equations. Proc. R. Soc. London, Ser. A 1969, 309, 195–208. (26) Avila, G.; Carrington, T., Jr. Non-Product Quadrature Grids to Solve the Vibrational Schroedinger Equation in 12D. J. Chem. Phys. 2009, 131, 174103. (27) Avila, G.; Carrington, T., Jr. Using Non-Product Quadrature Grids to Solve the Vibrational Schroedinger Equation in 12D. J. Chem. Phys. 2011, 134, 054126. (28) Kosloff, D.; Kosloff, R. A Fourier Method Solution for the Time Dependent Schroedinger Equation as a Tool in Molecular Dynamics. J. Comput. Phys. 1983, 52, 35–53. (29) Bramley, M. J.; Tromp, J. W.; Carrington, T., Jr.; Corey, G. C. Efficient Calculation of Highly Excited Vibrational Energy Levels of Floppy Molecules: The Band Origins of H+3 up to 35 000 cm1. J. Chem. Phys. 1994, 100, 6175–6194. (30) Chen, R.; Ma, G.; Guo, H. Six-Dimensional Quantum Calculations of Highly Excited Vibrational Energy Levels of Hydrogen Peroxide and its Deuterated Isotopomers. J. Chem. Phys. 2001, 114, 4763–4774. (31) Wang, X.-G.; Carrington, T., Jr. Six-Dimensional Variational Calculation of the Bending Energy Levels of HF Trimer and DF Trimer. J. Chem. Phys. 2001, 115, 9781–9796. (32) Wang, X.-G.; Carrington, T., Jr. Contracted Basis Lanczos Methods for Computing Numerically Exact Rovibrational Levels of Methane. J. Chem. Phys. 2004, 121, 2937–2954. (33) Chasman, D.; Beachy, M. D.; Wang, L.; Friesner, R. A. Parallel Pseudospectral Electronic Structure: I. Hartree-Fock Calculations. J. Comput. Chem. 1998, 19, 1017–1029. (34) Friesner, R. Solution of Self-Consistent Field Electronic Structure Equations by a Pseudospectral Method. Chem. Phys. Lett. 1985, 116, 39–43. (35) Ringnalda, M. N.; Won, Y.; Friesner, R. A. Pseudospectral Hartree-Fock Calculations on Glycine. J. Chem. Phys. 1990, 92, 1163– 1173. (36) Langlois, J.-M.; Muller, R. P.; Coley, T. R.; Goddard, W. A. I.; Ringnalda, M. N.; Won, Y.; Friesner, R. A. Pseudospectral Generalized Valence-Bond Calculations: Application to Methylene, Ethylene, and Silylene. J. Chem. Phys. 1990, 92, 7488–7497. (37) Friesner, R. Solution of the HartreeFock Equations for Polyatomic Molecules by a Pseudospectral Method. J. Chem. Phys. 1987, 86, 3522–3531. (38) Friesner, R.; Wyatt, R. E.; Hempel, C. A Generalized Version of the Recursive Residue Generation Method for Vector Computers. J. Comput. Phys. 1986, 64, 220–229. (39) Friesner, R. Solution of the HartreeFock Equations by a Pseudospectral Method: Application to Diatomic Molecules. J. Chem. Phys. 1986, 85, 1462–1468. (40) Martinez, T. J.; Carter, E. A. Pseudospectral Correlation Methods on Distributed Memory Parallel Architectures. Chem. Phys. Lett. 1995, 241, 490–496. (41) Beachy, M. D.; Chasman, D.; Friesner, R. A.; Murphy, R. B. Parallel Pseudospectral Electronic Structure: II. Localized MollerPlesset Calculations. J. Comput. Chem. 1998, 19, 1030–1038. (42) Murphy, R. B.; Beachy, M. D.; Friesner, R. A.; Ringnalda, M. N. Pseudospectral Localized Moller-Plesset Methods: Theory and Calculation of Conformational Energies. J. Chem. Phys. 1995, 103, 1481–1490. (43) Martinez, T. J.; Carter, E. A. Pseudospectral Multireference Single and Double Excitation Configuration Interaction. J. Chem. Phys. 1995, 102, 7564–7572. (44) Reynolds, G.; Martinez, T. J.; Carter, E. A. Local Weak Pairs Spectral and Pseudospectral Singles and Doubles Configuration Interaction. J. Chem. Phys. 1996, 105, 6455–6470. (45) Walter, D.; Szilva, A. B.; Niedfeldt, K.; Carter, E. A. Local WeakPairs Pseudospectral Multireference Configuration Interaction. J. Chem. Phys. 2002, 117, 1982–1993. (46) Kaminski, G. A.; Maple, J. R.; Murphy, R. B.; Braden, D. A.; Friesner, R. A. Pseudospectral Local Second-Order Moller-Plesset 2198

dx.doi.org/10.1021/jz200513h |J. Phys. Chem. Lett. 2011, 2, 2193–2199

The Journal of Physical Chemistry Letters Methods for Computation of Hydrogen Bonding Energies of Molecular Pairs. J. Chem. Theory Comput. 2005, 1, 248–254. (47) Greeley, B. H.; Russo, T. V.; Mainz, D. T.; Friesner, R. A.; Langlois, J.-M.; Goddard, W. A. I.; Donnelly, R. E.; Ringnalda, M. N. New Pseudospectral Algorithms for Electronic Structure Calculations: Length Scale Separation and Analytic Two-Electron Integral Corrections. J. Chem. Phys. 1994, 101, 4028–4041. (48) Murphy, R. B.; Cao, Y.; Beachy, M. D.; Ringnalda, M. N.; Friesner, R. A. Efficient Pseudospectral Methods for Density Functional Calculations. J. Chem. Phys. 2000, 112, 10131–10140. (49) Ko, C.; Malick, D. K.; Braden, D. A.; Friesner, R. A.; Martinez, T. J. Pseudospectral Time-Dependent Density Functional Theory. J. Chem. Phys. 2008, 128, 104103. (50) Becke, A. D.; Dickson, R. M. Numerical Solution of Shrodinger’s Equation in Polyatomic Molecules. J. Chem. Phys. 1990, 92, 3610–3612. (51) Dickson, R. M.; Becke, A. D. Basis-Set-Free Local DensityFunctional Calculations of Geometries of Polyatomic Molecules. J. Chem. Phys. 1993, 99, 3898–3905. (52) Chelikowsky, J. R.; Troullier, N.; Wu, K.; Saad, Y. HigherOrder Finite-Difference Pseudopotential Method - An Application to Diatomic-Molecules. Phys. Rev. B 1994, 50, 11355–11364. (53) Briggs, E. L.; Sullivan, D. J.; Bernholc, J. Real-Space MultigridBased Approach to Large-Scale Electronic Structure Calculations. Phys. Rev. B 1996, 54, 14362–14375. (54) Castro, A.; Appel, H.; Oliveira, M.; Rozzi, C. A.; Andrade, X.; Lorenzen, F.; Marques, M. A. L.; Gross, E. K. U.; Rubio, A. OCTOPUS: A Tool for the Application of Time-Dependent Density Functional Theory. Phys. Stat. Sol. B 2006, 243, 2465–2488. (55) Marques, M. A. L.; Castro, A.; Bertsch, G. F.; Rubio, A. OCTOPUS: A First-Principle Tool for Excited ElectronIon Dynamics. Comput. Phys. Commun. 2003, 151, 60–78. (56) Grabowski, P. E.; Chernoff, D. E. Pseudospectral Calculation of the Wave Function of Helium and the Negative Hydrogen Ion. Phys. Rev. A 2010, 81, 032508. (57) Hassoun, M. H. Fundamentals of Artificial Neural Networks; MIT Press: Cambridge, MA, 1995. (58) Lagaris, I. E.; Likas, A.; Fotiadis, D. I. Artificial Neural Network Methods in Quantum Mechanics. Comput. Phys. Commun. 1997, 104, 1–14. (59) Sugawara, M. Numerical Solution of the Schrodinger Equation by Neural Network and Genetic Algorithm. Comput. Phys. Commun. 2001, 140, 366–380. (60) Nakanishi, H.; Sugawara, M. Numerical Solution of the Schr€odinger Equation by a Microgenetic Algorithm. Chem. Phys. Lett. 2000, 327, 429–438. (61) Manzhos, S.; Carrington, T., Jr. An Improved Neural Network Method for Solving the Schrodinger Equation. Can. J. Chem. 2009, 87, 864–871. (62) Caetano, C.; Reis, J. L., Jr.; Amorim, J.; Ruv Lemes, M.; Dal Pino, A., Jr. Using Neural Networks to Solve Nonlinear Differential Equations in Atomic and Molecular Physics. Int. J. Quantum Chem. 2010, 111, 2732–2740. (63) Manzhos, S.; Yamashita, K.; Carrington, T., Jr. Calculating Anharmonic Vibrational Frequencies of Molecules Adsorbed on Surfaces Directly from Ab Initio Energies with a Molecule-Independent Method: H2O on Pt(111). Surf. Sci. 2011, 605, 616–622. (64) Manzhos, S.; Yamashita, K.; Carrington, T. Using a Neural Network Based Method to Solve the Vibrational Schr€odinger Equation for H2O. Chem. Phys. Lett. 2009, 474, 217–221. (65) Manzhos, S.; Carrington Jr., T.; Yamashita, K. Solving the Vibrational Schr€odinger Equation without a Potential Energy Surface Using a Combined Neural Network Collocation Approach. In Computer Physics; Doherty, B. S., Molloy, A. N., Eds.; Nova Science Publishers, Inc.: Hauppauge, NY, 2011. (66) Nakatsuji, H.; Nakashima, H.; Kurokawa, Y.; Ishikawa, A. Solving the Schr€odinger Equation of Atoms and Molecules without Analytical Integration Based on the Free Iterative-Complement-Interaction Wave Function. Phys. Rev. Lett. 2007, 99, 240402.

PERSPECTIVE

(67) Bande, A.; Nakashima, H.; Nakatsuji, H. LiH Potential Energy Curves for Ground and Excited States with the Free Compliment Local Schr€odinger Equation Method. Chem. Phys. Lett. 2010, 496, 347–350. (68) Bramley, M. J.; Carrington, T., Jr. A General Discrete Variable Method to Calculate Vibrational Energy Levels of Three- and FourAtom Molecules. J. Chem. Phys. 1993, 99, 8519–8541. (69) Wang, X.-G.; Carrington, T., Jr. Parallel Methods for HighDimensional Quantum Dynamics. Comput. Phys. Commun. 2010, 181, 455–461. (70) Parallel Computing in Chemical Physics, 128th ed. Comput. Phys. Commun. 2000, 1530.

2199

dx.doi.org/10.1021/jz200513h |J. Phys. Chem. Lett. 2011, 2, 2193–2199