Article pubs.acs.org/jchemeduc
A Systematic Approach for Understanding Slater−Gaussian Functions in Computational Chemistry Brianna Stewart, Derrick J. Hylton, and Natarajan Ravi* Physics Department, Spelman College, Atlanta, Georgia 30314, United States S Supporting Information *
ABSTRACT: A systematic way to understand the intricacies of quantum mechanical computations done by a software package known as “Gaussian” is undertaken via an undergraduate research project. These computations involve the evaluation of key parameters in a fitting procedure to express a Slater-type orbital (STO) function in terms of the linear combination of Gaussian-type orbital (GTO) functions. A procedure for the optimization process based on the Newton−Raphson method is developed and is applied to STO-2G and STO-3G basis sets. Satisfactory results obtained by this procedure are used to illustrate the importance of ab initio computations for inclusion in the chemistry or physics undergraduate curriculum. Programming languages such as Python and Maple were employed to obtain the results. KEYWORDS: Upper-Division Undergraduate, Graduate Education/Research, Physical Chemistry, Computer-Based Learning, Computational Chemistry, Mathematics/Symbolic Mathematics, Quantum Chemistry
A
will lead to a better understanding of one of the most used quantum-mechanical software package “Gaussian”. Toward achieving this goal, the student was trained to perform all the necessary numerical integrals using Python programming and cross-check the results with Maple. Hence, we set out to investigate this problem and developed a new procedure to evaluate α with goals of (1) helping students understand the intricacies of the computation and (2) disseminating the procedure to the community of the Gaussian program users. The developed procedure employs the Newton−Raphson method for multivariable parameters to achieve the minimization process. In this manuscript, the theory, an example of the method, results, and applications are presented.
majority of undergraduate students, upon completion of the first level of a quantum mechanics course, are unsure of the concepts learned and grudge the mathematics involved in understanding the course material. The original lack of understanding of the fundamentals prevents students from moving forward to develop an algorithm for a given situation. In spite of the mathematical drudgery, motivated students, when undertaking research in the field of quantum computation, find it hard to comprehend the level of math involved, but only gradually develop an appreciation for the amount of work needed to address a research question. In modern times, it is true that a sizable number of physics and chemistry students conducting undergraduate research in computational fields use commercial software package ranging from worksheet program such as Excel to quantum mechanical programs such as Gaussian.1 Using such built-in programs offers flexibility and ease, on the one hand but, on the other hand, obscures the mathematical details that take place behind the scene. To bridge the gap between conceptual understanding and using commercial programs as black boxes, we set out to perform computations in the conventional sense as a first step that involves evaluating necessary integrals, writing codes for appropriate situations, and making sense of the results as an undergraduate research project. In particular, this project involved computing the coefficients that are already built-in in the Gaussian software package. A detailed procedure to evaluate the key exponential parameter α and the coefficients that connect a Slater-type function with a Gaussian-type function was not available in the literature or, when available, the details were sketchy at best.2−4 This attempt, we believe, © XXXX American Chemical Society and Division of Chemical Education, Inc.
■
THEORY The application of quantum mechanical methods to investigate molecular properties involves sophisticated calculations and requires a detailed algorithm or program to translate the mathematical formalism into “step-by-step” instructions suitable for execution on a computer. Gaussian is a commercial software package that is used today to achieve this objective. Complex molecular structures are studied to yield a wealth of information by using this approach. For example, calculations using Gaussian for Fe−sugar complex system have been reported in predicting the geometries and energies and are correlated with experimental data. Similar calculations on methane intercalated graphene layers offer condition to
A
dx.doi.org/10.1021/ed300807y | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
Article
experimentally prepare a tunable graphene based semiconductor.5 The fundamental principle in this method is to express a Slater-type orbital (STO) in terms of Gaussian functions followed by an optimization procedure. Normalized GTO STO, φSTO 1s , and Gaussian-type orbital (GTO) functions, φ1s , are expressed as follows for a 1s orbital ⎛ ζ 3/2 ⎞ ⎟e−ζr ⎜ φ1STO = s ⎝ π ⎠ φ1GTO = s
We note that if the Slater orbital is not an s-type, evaluating the angular part of the integral is not as straightforward as it is for an s-type orbital.
■
METHOD FOR STO-2G To illustrate the method, we consider STO-2G, wherein there will be two Gaussian terms with two different α values (α1 and α2) and two coefficients (C1 and C2) corresponding to each Gaussian term. We need to determine these four parameters, subject to the normalization constraint, to minimize eq 5 or maximize the sum of terms like eq 6. The first step is to use the normalization constraint to write the second integral of eq 5 as unity and to eliminate one of the parameters, so that we are left with three independent parameters. With this idea, eq 5 for N = 2 becomes
(1)
⎛ 2α ⎞3/4 −αr 2 ⎜ ⎟ e ⎝π ⎠
(2)
where r is the spherical coordinate and α and ζ are adjustable parameters; for s-type orbitals, there is no angular dependence. The idea that integrals involving Slater-type functions are timeconsuming and rather difficult to evaluate, due to slow convergence, led scientists in the 1950s to express Slater-type functions in terms of Gaussian functions as6 ⎛ ζ 3/2 ⎞ ⎟ e−ζr = ⎜ ⎝ π ⎠
N
∑ i=0
2 ⎛ 2α ⎞3/4 Ci⎜ i ⎟ e−αir ≡ φ1NG s ⎝ π ⎠
∫
=2−2
(3)
(φ1STO − φ1NG )2 dτ s s
1=
∫
dτ +
∫
dτ − 2
∫
φ1STO φ1NG s s
∫0
∞
(φ1GTO (α1))2 dτ s
φ1GTO (α1)φ1GTO (α 2 ) d τ s s
(φ1STO (α2))2 dτ = C12 + 2C1C2K (α1 , α2) s
K(α1 , α2) ≡
(9)
∫
φ1GTO (α1)φ1GTO (α 2 ) d τ = s s
(4α1α2) 3/4 (α1 + α2) 3/2 (9a)
Equation 9 is treated as a quadratic equation in C2 and solved to get (positive solution) C2 =
dτ
C12K 2 − C12 + 1 − C1K
(10)
where K = K(α1,α2). With the elimination of this parameter, we write the square brackets of eq 8 as the quantity to be maximized in order to minimize the least-squares integral F(α1 , α2 , C1) ≡ C1S(α1) + ( C12K 2 − C12 + 1 − C1K ) S (α 2 )
(11)
To achieve the optimization, one imposes the optimization condition by differentiating the expressions with respect to each parameter. Equating each differentiated expression separately to zero to achieve the optimization condition, the unknown parameters can be determined using the Newton−Raphson (N−R) method for three variables. Instead of the N−R method, one can use Microsoft Excel Solver for the optimization process.9 However, our emphasis is for students to understand the underlying mathematics, and hence, we refrained from using a black box approach and resorted to the
(6)
2
r 2e−ζr e−αir dr
∫
∫
∫
where
When one uses spherical polar coordinates, and taking into account the angular part that is evaluated to be 4π for 1s orbitals, the integral that consists solely of the radial part is ⎛ ζ ⎞3/2 ⎛ 2α ⎞3/4 S(αi) = 4π ⎜ ⎟ ⎜ i ⎟ ⎝π ⎠ ⎝ π ⎠
(8)
+ C22
where the first term is evaluated to be unity because φSTO is 1s normalized, and the third integral is known as the overlap integral, which should be maximum for the optimization because of the negative sign so that eq 5 becomes minimum. The overlap integral can be evaluated numerically and consists of the sum of N terms of the form
∫
(φ12sG)2 dτ = C12
+ C22
(5)
⎛ ζ 3/2 ⎞ 2 ⎛ 2α ⎞3/4 ⎟ e−ζr ⎜ i ⎟ e−αir dτ ≡ S(αi) ⎜ ⎝ π ⎠ ⎝ π ⎠
∫
+ 2C1C2
(4)
φ1NG φ1NG s s
φ1STO φ12sG dτ s
where S(α) is defined in eq 7. We use the normalization constraint for the elimination of one of the parameters, C2, as follows:
which is the least-squares fitting criterion for a function, where dτ represents the differential volume element. The minimization is done subject to the constraint that φNG 1s is normalized. Expansion of this integral results in three integrals, φ1STO φ1STO s s
∫
= 2 − 2[C1S(α1) + C2S(α2)]
where N ≥ 2 is a chosen integer that determines how many GTO will be used. In other words, a linear combination of Gaussian functions was made to be the appropriate eigenfunctions of a Hamiltonian for which Slater-type functions may be exact eigenfunctions, as suggested by the hydrogen problem. The important parameters are αi and Ci that are to be determined to satisfy eq 3. The convention STO-3G, for instance, which forms one of the basis sets in the Gaussian program, denotes that a Slater function is expressed as a sum of three Gaussian terms with varying α and C values. Usually, the procedure to extract these α values is done by an optimization process of minimizing the following integral:7,8
∫
− φ1NG (φ1STO )2 dτ s s
(7) B
dx.doi.org/10.1021/ed300807y | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
Article
functions. For instance, an attempt to find the parameters for a STO-3G scenario resulted in values that are in agreement with the reported values. However, the procedure also converged to 1G and 2G solutions, as these are also possible solutions for the 3G problem. This is due to the fact that there are as many minima as the number of Gaussian terms; reaching the third narrow minimum is always difficult and mathematics facilitates reaching more easily the broader minimum rather than the narrower minimum. To arrive at a solution, the starting parameters must be somewhat closer to the real 3G values. Based on an educated guess or by trial and error, it is not difficult to start with reasonably correct values for the iterative process. For instance, for all α values larger than 0.2, a nonconverging solution was seen in the 3G problem, which amounted to one of the α values being less than 0.2. For N > 3, because of the difficulty of the method, other local or global optimization methods may be more effective.11,12 This basic idea of expanding a STO in terms of a Gaussian basis set plays a central role in not only Gaussian software package, but also in more sophisticated computations using unrestricted Hartree−Fock theory (UHF) and density functional theory (DFT).13−15 It is also noteworthy to mention that one can follow a procedure based on Lagrange multiplier, which also gives rise to the same values. The salient features of this, although described generally, are available in the literature in a paper published in 1969.4,9 For the benefit of the readers, we have reproduced the necessary mathematical derivations and other details in the Supporting Information.
N−R method. The N−R method is used to solve the system of equations ∂F ∂F ∂F = =0 = ∂α1 ∂α2 ∂C1
(12)
Defining ∂F/∂α = f1, ∂F/∂α2 = f 2, ∂F/∂C1 = fc, the N−R method is a procedure to correct the previous approximate solution in an iterative process.10 The correction is determined by −1 ⎛ ∂f ∂f1 ∂f1 ⎞ 1 ⎟ ⎜ ⎜ ∂α ∂α2 ∂c1 ⎟ ⎛ − f (α , α , C ) ⎞ ⎛ Δα1 ⎞ ⎜ 1 2 1 1 1 ⎟ ⎟ ⎟ ⎜ ∂f ⎜ ∂f2 ∂f2 ⎟ ⎜⎜ 2 ⎟ f C α α ( , , ) − ⎜ Δα2 ⎟ = ⎜ 2 1 2 1 ⎟ ⎟ ⎟ ⎜ ∂α1 ∂α2 ∂c1 ⎟ ⎜ ⎜ ⎜− f (α , α , C )⎟ ⎝ ΔC1 ⎠ ⎜ 1 2 1 ⎝ ⎠ C ∂f ∂fc ∂fc ⎟ ⎟⎟ ⎜⎜ c ⎝ ∂α1 ∂α2 ∂c1 ⎠
(13)
where Δα1, Δα2, ΔC1 are the corrections to the previous values of the approximate solutions to eq 12. To get expressions for these matrix elements, one differentiates the function in eq 11 to get ∂K ⎛ ⎞ Kc12 ∂α dS ∂K ⎟ ⎜ 1 − c1 f1 = c1 (α1) + ⎜ S (α 2 ) = 0 2 2 2 ∂α1 ⎟ dα − + 1 c K c 1 ⎝ 1 ⎠
(14)
dS (α 2 ) dα ∂K ⎞ Kc12 ∂α ∂K ⎟ 2 − c1 S( α 2 ) = 0 ∂α2 ⎟ (c12K 2 − c12 + 1) ⎠
f2 = ( (c12K 2 − c12 + 1) − c1K ) ⎛ ⎜ +⎜ ⎝
⎡ ⎤ C1K 2 − C1 ⎢ fc = S(α1) + − K ⎥S(α 2 ) = 0 ⎢ C 2K 2 − C 2 + 1 ⎥ ⎣ ⎦ 1 1
■
(15)
(16)
If Maple is used, the integration in S(α) and dS/dα are written as an integral that are evaluated internally; if Python is used, a numerical integration program needs to be written. It is straightforward to calculate ∂K/∂α1 and ∂K/∂α2 from eq 9a. We derive an expression for (dS/dα) from eq 7 to yield ⎛ζ ⎞ dS = 4π ⎜ ⎟ ⎝π ⎠ dα 0
−
3/2
⎛ 2α ⎞3/4 ⎡ 3 ⎟ ⎝ π ⎠ ⎢⎣ 4α
⎜
∫∞ r 4e−ζre−α r
2
0
∫∞ r 2e−ζre−α r
2
dr
⎤ dr ⎥ ⎦
(17)
For the derivatives in eq 13, we can use an approximation such as ∂f f (α1 + ε) − f (α1) = ∂α1 ε
where ε = 10
■
−7
APPLICATIONS
Building computational skills in chemistry curriculum is extremely important and is currently addressed in chemistry education research. We suggest introducing exercises involving computations in chemistry from general chemistry onward and using Maple may be an effective tool to achieve this objective. The article illustrates this idea. Although the current problem is not for general chemistry, Maple is an effective tool to build computational skills among students. The method enumerated above may be used for various situations encountered in chemistry or chemical education. We list a few plausible scenarios: (a) Computations that require optimization procedure can use the Newton−Raphson method described in this article or Lagrange multiplier method given in the Supporting Information and apply it to this problem as a nontrivial application. (b) These computational techniques can be taught in several courses such as computational chemistry or physical chemistry in a way that connects with modern topics in chemistry. (c) This method demonstrates one potential use of Maple in teaching chemistry, and we recommend programming and algorithm development as mechanisms to enhance computational skills. This method is particularly suitable for undergraduate research projects and some suggested undergraduate research projects are (i) extend this approach using the 2s orbital to obtain the coefficients of the Gaussian expansion, (ii) determine the equilibrium concentration of species given the equilibrium constant and the initial concentrations, and (iii) determine the first few coefficients of the Fourier expansion via this technique presented here to expand a given function in terms of a set of functions.
(18)
for a four figure accuracy on the solution.
RESULTS This procedure yielded α and C values in agreement with the reported values.2 However, it should be pointed out that, although this procedure gives satisfactory results, it poses more challenges for situations that involve more than two Gaussian C
dx.doi.org/10.1021/ed300807y | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
■
Article
Publishers: New York, 2011; pp 1−48. (b) Hargrove, J.; Mihiri Shashikala, H. B.; Guerrido, L.; Ravi, N.; Wang, X.-Q. Band Gap Opening in Methane Intercalated Graphene. Nanoscale 2012, 4, 4443−4446. (6) Boys, S. F. Electronic Wave Functions I. A General Method of Calculation for the Stationary States of Any Molecular System. Proc. R. Soc. London 1950, 542−554. (7) Hehre, W. J.; Radom, L.; Schleyer, P.v.R.; Pople, J. A. Ab Initio Molecular Orbital Theory; John Wiley & Sons: New York, 1986; p 68. (8) Lewars, E. Computational Chemistry Introduction to the Theory and Applications of Molecular and Quantum Mechanics, Kluwer Academic Publishers (2003) p188. (9) Pye, C. C.; Mercer, C. J. On the Least-Squares Fitting of Slatertype Orbitals with Gaussians: Reproduction of the STO-NG Fits using Microsoft Excel and Maple. J. Chem. Educ. 2012, 89 (11), 1405−1410. (10) (a) Module for Fixed Point Iteration and Newton’s Method in 2D and 3D. http://math.fullerton.edu/mathews/n2003/ FixPointNewtonMod.html and Internet resources and bibliography contained therein (accessed Mar 2013). (b) Chan, T. F. An Approximate Newton Method for Coupled Nonlinear Systems. SIAM J. Numer Anal. 1985, 22 (5), 904−913. (11) Handbook of Global Optimization; Pardalous, P. M., Romeijn, H. E. , Eds.; Kluwer Academic Publishers: Boston, MA, 2002; Vol. 2. (12) Schlick, T. Optimization Methods in Computational Chemistry. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; John Wiley and Sons, Inc.; Hoboken, NJ, 2007; Vol. 3. (13) (a) Becke, D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. 1988, A38, 3098−3100. (b) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. 1988, B37, 785−789. (14) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. 12. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular-Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257−2262. (15) Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-Consistent Molecular Orbital Methods. 9. Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules. J. Chem. Phys. 1971, 54, 724−728.
CONCLUSIONS The following main conclusions are derived from the present study. 1. An optimization procedure using the Newton−Raphson method for multivariable parameters can be used to evaluate the key parameters α and the coefficients for a Gaussian expansion of a Slater orbital, after using the normalization constraint to eliminate one of the coefficients. 2. This procedure can also be extended to 4−6 G basis functions for further computations, although the minima are narrow and may not be easy to detect. 3. This procedure can be a useful tool to be introduced to students and may constitute an important curricular exercise for computational undergraduate research.
■
ASSOCIATED CONTENT
S Supporting Information *
Lagrange multiplier method. This material is available via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS N.R. acknowledges the support by the National Science Foundation Grant PREM DMR-093412 and the Department of Energy Grant DE-FG52-09NA29518. The authors would like to thank Beatriz Cardelino for useful discussions and helpful suggestions.
■
REFERENCES
(1) 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., Nakatsuji, H., Caricato, M., Li, X., Hratchian, H. P., Izmaylov, A. F., Bloino, J., Zheng, G., Sonnenberg, J. L., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Vreven, T., Montgomery, Jr., J. A., Peralta, J. E., Ogliaro, F., Bearpark, M., Heyd, J. J., Brothers, E., Kudin, K. N., Staroverov, V. N., Kobayashi, R., Normand, J., Raghavachari, K., Rendell, A., Burant, J. C. Iyengar, S. S. Tomasi, J. Cossi, M. R., Millam, N. J., Klene, M. Knox, J. E., Cross, J. B., Bakken, V., Adamo, C., Jaramillo, J., Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. Martin, R. L., Morokuma, K., Zakrzewski, V. G., Voth, G. A., Salvador, P., Dannenberg, J. J., Dapprich, S., Daniels, A. D., Farkas, O., Foresman, J. B., Ortiz, J. V., Cioslowski, J., Fox, D. J. Gaussian 09, Revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (2) Hahre, W. J.; Stewart, R. F.; Pople, J. A. Self Consistent Molecular Orbital Methods. I. Use of Gaussian Expansions of Slater Type Atomic Orbitals. J. Chem. Phys. 1969, 51, 2657−2664. (3) Hahre, W. J.; Ditchfield, R.; Stewart, R. F.; Pople, J .A. SelfConsistent Molecular Orbital Methods. IV. Use of Gaussian Expansions of Slater Type of Orbitals. Extentsion to Second Row Molecules. J. Chem. Phys. 1970, 52, 2769−2773. (4) Stewart, R. F. Small Gaussian Expansions of Atomic Orbitals. J. Chem. Phys. 1969, 50, 2485−2495. (5) (a) Jeilani, Y.; Cardelino, B. H.; Ravi, N. 57Fe Mössbauer Spectroscopic and Density Functional Theory (DFT) Study on the Interactions of the Metal Ion with Monosaccharides; Nova Science D
dx.doi.org/10.1021/ed300807y | J. Chem. Educ. XXXX, XXX, XXX−XXX
