Computers Chem. Vol. I?, No. 3, pp. 287-290, 1993 Printed in Great Britain. All rights reserved
0097.8485/93 56.00+ 0.00 Copyright 0 I993 Pergamon Press Ltd
CALCULATION OF TWO-CENTER INTEGRALS BETWEEN SLATER-TYPE ORBITALS
R. WOJNECKI’* and P. MODRAK~ ‘Institute of Physics of the Polish Academy of Sciences, Al. Lotnik6w 32/46, 02-668 Warszawa and
*Institute of Physical Chemistry of the Polish Academy of Sciences, Kasprzaka 44/52,01-224 Warszawa, Poland
(Received 6 May 1992; in revised form 9 February 1993)
Abstract-Two-center overlap, nuclear-attraction and kinetic energy integrals between Slater-type orbitals (STOs) centered on two different points are calculated. The expansion of ST0 in terms of spherical harmonics centered on an atom dispiaoed from the first one is used.~CoelTicients of the expansion are stored in an adjustable array (its size depends on quantum numbers of the ST0 involved) in the first step of the calculation. In the second step a simple and closed formula is used for obtaining the required type of integral between STOs (with quantum numbers in the range for which calculations in the first step have been performed).
INTRODUCTION
Slater type orbit& (STOs) constitute a very con- venient set of orbitals for calculations of many elcc- tronic properties of molecules and solids (Davenport, 1984; Fernandez Rico et al., 1988; Guseinov, 1988; Jones, 1988; Kotani et al., 1955; Modrak & Woj- necki, 1987, 1989; Roothaan, 1951, 1956; Rueden- berg, 1951; Ruedenberg et al., 1956; Sharma, 1976; SufXczyAski, 1956; Weinert et al., 1986; Watson ef al., 1986, 1987). The calculation of the LGwdin a function and the integrals between STOs has been recently subject of many papers (Jones, 1992; Homeier & Steinbom, 1992 and references therein). The purpose of the present paper is to present the efficient program for calculating integrals between STOs which can be useful for many chemists and physicists.
The program described in the present paper is able to compute two-center overlap, nuclear attraction and kinetic energy integrals between STOs for any set of quantum numbers. The program is designed to be a subroutine of a main program. In this form it is easy to use in any kind of calculations where integrals between STOs are required.
The formulas for calculations of the integrals are taken from Sharma (1976). They were checked by the authors and limits of summations in equation (23~) of Sharma (1976) were corrected [to avoid (-n)! in the formula]. The sum in equation (36) of Sharma (1976) was rearranged to avoid serious numerical errors as described in the comments enclosed in the description of the function GN.
In order to test the calculations, the results of the present program have been compared with the results obtained with the use of the formulas given by
* Author for correspondence.
Roothaan (1951) and Suffczyfiski (1956). The agree- ment is excellent, relative error is < lo-“. The values of the test parameters were taken close to the values used in practical calculations (Modrak et al., 1987). One can also teat the accuracy of the program for any set of the parameters with the aid of the supplied test program.
DESCRIPI-ION OF THE PROGRAM
The program consists of two main subroutines and six functions. The segment playing the role of the main program (DEMO), the test program (TESTSL) and two subroutines for test calculations can be also supplied. The organization of the program with the use of test programs as a main programs is shown in Fig. 1. In the figure NAME (written in upper-case letters) is the name of the program, subroutine or function, name (written in lower-case letters) is the name of the file produced by the program. The segments playing the role of the main program is either TESTSL or DEMO. TAB and SLINT are two main subroutines, BVSLLM, FKK, BET, CAN, GN, HN are functions. The program TESTSL makes use of two additional subroutines, ROOT and SUFF, to test results and creates the file testsims to store results of the test.
Integrals are calculated for the product of the two following functions:
fn = ANORM rf”exp(AM * rA) * YLM.MC(v, cp)
(1)
and
fa = ANORD riDexp(AD * r,) 8 Y,,.,,(tl, 4)
(2)
287
288 R. Wo~mcn and P. MODRAK
where r,, v, p are the spherical coordinates expressed in the frame of reference centered at the point A, rs, 8, q5 are the spherical coordinates expressed in the frame of reference centered at the point B; FL,, are the spherical harmonics; NM, AM, LM, M and ND, AD, LD, MC are the program variables defining the function form off., and fs , respectively (both AM and AD must be negative). The norms ANORM and ANORD are calculated in the program. Tbe point B has (0, 0, DZST) coordinates in the Cartesian axes system centered at the point A (DZST is the program variable). The function fB can be expressed in terms of the spherical harmonics centered at the point A using formulae given in Sharma (1976). The overlap integral between unnormalized STOs takes finally the following form (Sharma, 1976):
(NM, LM, MC, AM, 1 ND, LD, MC, AD) =
GN(N, X) = N!/XN+’ - exp( - X) * iYN(N, X) (7)
FKI.x(ND, LM, LD, MC) SMAX VMAX
= &vt v-c,, BV(S, LM, LD, MC)
42 + S - LD + ND)!/[(KZ - 2 * V)!
*(2*S-LD+ND-K-KI+2*V)!] (8)
SMIN = MAX(0, {l/2 * (K + LD - ND + 1))) (9)
VMIN = MAX(0, {l/2 * (KI + K
+LD-ND-2*S+l))) (10)
VMAX = MIN(LM + LD - S, {l/2 * KI}).
(11)
{a} is the integer part of a, f3V (S, LM, LD, M) is given by equation (17b) in Sharma (1976).
DZSTND + NM + k Nlr [exp(AD * DZST)
x (-NS! * [DZST L (AM + AD)]-Ns- ’
+ (- l)U * GN[NS, DIST * (AD - AM)]) KMAX
x c (-AD mDZST)-K-’ K-=0
* FKI.x(ND, LM, LD, MC)
-exp(AM + DZST) + (- I)LDf“‘D-u
* HN[NS, DZST * (-AM - AD)]
KMAY x c (AD * DZST)+’
KC0
* F&ND, LM, LD, MC) 1 where:
NS=NM-LM+KZ
KMAX=2*LM+LD+ND_KI
The main program
The main program should contain:
(i) The statement PARAMETER defining values of the following variables: LMMIN, LMMAX, LD- MIN, LDMAX, NDMZN, NDMAX, LFI, LF2, NSZL, MM and NFKT. LMMIN and LMMAX give the range of LM [see equation (l)] in which the calculations are to be performed. LDMZN and LD- MAX define the range of LD [see equation (2)]. NDMZN and NDMAX define the range of ND [see equation (2)]. The parameters defined above and the maximum value NMMAX of NM [see equation (1)] which is to be used in the intended calculation
(3) constitute the set of the independent parameters which are fixed by the user. The remaining par- ameters depend on this set in the following way:
(4) MM = MIN(LDMAX, LMMAX)
(5) LFl = LMMAX + LDMAX
LF2 = NDMAX + 2 * LMMAX + LDMAX (6)
NSIL = MAX(LF2, NMMAX + LDMAX).
TEST.9 DEMO I
I ROOT SUFF TAB
test sl . res
FKK
I ’
CAN I
I I BVSLLM GN !-IN
Fig. 1. The organization of the program.
Two-center integrals between SLOs 289
NFKT is the number of the non-zero coefficients FKIX (ND, LM, LD, MC) [see equation (811 for the given range of quantum numbers. NFKT =820 for the range used in the program DEMO and TESTSL, that is for:
LMMIN = 0; LMMAX = 2; LDMIN = 0;
LDMAX = 2; NDMIN = 1; NDMAX = 3;
NMMAX = 3.
Any other program can use the same value of NFKT provided that the range of quantum numbers does not exceed the range given above. The recipe to get the proper value of NFKT (minimizing a computer memory required) is the following:
(a) put NFKT = 1,
(b) call the subroutine TAB.
The subroutine wil1 write the proper value of NFKT and stop afterwards.
(ii) The dimension statements:
INTEGER * 2 ZFK(O:LFl, O:LF2, NDMIN:NDMAX, LMMIN:LMMAX, LDMIN:LDMAX, &MM) REAL * 8 SIL(O:NSIL), FKT(O:NFKT).
(iii) The call statements:
CALL TAB (LMMIN, LMMAX, LDMIN, LDMAX, NDMIN, NDMAX, NMMAX, IFK, EKT, NFKT, SIL, NSIL, LFl, LF2, MM).
The subroutine TAB needs to be called only once for all integrals involving STOs with the parameters LM, NM and LD, ND confined to the range defined by the first seven parameters of the subroutine TAB.
CALL SLINT (IND, NM, ND, LM, LD, AM, AD, DIST, MC, RES, LMMIN, LMMAX, LDMIN, LDMAX, NDMIN, NDMAX, NMMAX, IFK, FKT, NFKT, SIL, NSIL, LFl, LF2, MM)
where:
IND is the parameter describing the type of the calculated integral: for IND = 1 the overlap integral is calculated; for IND = 2 the nuclear-attraction inte- gral is calculated; for IND = 3 the kinetic energy integral is calculated.
NM, ND, LM, LD, AM, AD, MC describe the Slater orbitals involved [see equations (1) and (2)], DIST is the distant between displaced centers and the output value of RES gives the desired value of the integral.
The program DEMO gives the value of a chosen integral for a given set of parameters. The program TESTSL also compares integrals calculated by the present method for a given set of parameters with the integrals calculated with the aid of formulae given by Roothaan (1951) and Suffczytiski (1956) in all cases for which the latter formulae are applicable, writes
results of this comparison to the file (tests&s) and the result with the maximum error to the screen.
The relative error is usually very small, but grows with the decrease of the value of DIST a (-AM) and DIST * (-AD). The message “IMPROPER PARAMETERS” is given and the program stops when either of these values is qO.1.
Subroutine TAB
The meaning of the parameters of the subroutine TAB has been described in the preceding section.
Description of the subroutine. In this subroutine ali parameters depending on a chosen range of quantum numbers describing adjustable array dimensions are checked. In the case of an error the message is written to the screen:
“parameter name” = “proper value” NOT “present value”
and the program stops. If all parameters are correct the values of n! are calculated and stored in the SIL array. Then, the coefficients FKIK (ND, LM, LD, M) given by equation (8) in the range of KI and K defined by the limits of sums in equations (3) and (5) are calculated with the aid of the function FKK. To save computer memory only nonzero coefficients are stored in the array FKT. The value of IFK(KI, K, ND, LM, LD, M) gives the index of the matrix element FKT storing value of the coefficient F,, (ND, LM, LD, M) (the index zero corresponds to the coefficients equal to zero).
The function FKK makes the use of the function BVSLLM to calculate the coefficients b,(sfLM) given by equation (17b) in Sharma (1976).
Subroutine SLINT
The meaning of the parameters of the subroutine SLINT has been already described in Main Program.
Description of the subroutine. First, the correctness of parameters is checked. If either of DIST, DIST * (-AM), or DIST * (- AD) is too small; AM or AD is ~0; either of the parameters LM, LD, NM, ND is out of the range; ABS (MC) is greater than the minimum of LM and LD values, a message “IM- PROPER PARAMETERS” is written to the Screen and the program stops. If ah parameters have the proper values, the norms of the STOs are calculated and the “IF” statement makes the program to calcu- late the desired type of the integral defined by the value of the parameter IND:
(a) Overlap integrals (IND = 1). The result is equal to the product of norms of STOs multiplied by the overlap integral between two unnormalixed STOs [given by equations (1) and (2)] calculated with the aid of the function CAN.
(b) Nuclear-attraction integrals (IND = 2). The function from equation (1) is multiplied by (-l/r,). One can write the result as:
fA 3 -,iM-’ ew(AM * rA) * YLM.MC(~I cp>. (12)
290 R. WOJNECKI and P. MODRAK
Then the function CAN is used with NM - 1 par- ameter in the place of NM and the result is multiplied by norms of STOs and taken with minus sign.
(c) Kinetic energy integrals (fND = 3). The kinetic energy operator - l/2 A in spherical coordinates ap plied to the function given by equation (1) gives:
-1/2A(f,)=[-l/2cAM2*rzM-2*AAM*NM
*r~“-‘-(NM*(Nh4-l)-LM
* (LM + 1)) * rh-‘1
*e&AM + rA) * YL.+f.yc(v, cp). (13)
The function CAN is used with the parameter NM, then with NM - 1, and finally with NM - 2 and multiplied by the norms and the proper factors, appearing in equation (13).
The function CAN makes use of the function B VSL LM.
Function BVSLLM
The subroutine calculates values b,(sZLM), making use of the function BET to calculate the factor following the square root in equation (17~) of Sharma (1976).
Function BET
Calculates the factor following the square root in equation (17~) of Sharma (1976).
Function FKK
Calculates the expression occuring in equation (8) with corrected lower limits of the summations [see equations (9), (lo)].
Function CAN
Calculates the integral between two STOs given by equations (1) and (2) using formula (3). Two func- tions GN and IfN are called during calculations and the arrays IFK and FKT are used.
Function GN
Calculates the value of GN (N, X) given by equation (7). The equation is equivalent to equation (35d) of Sharma (1976). This equation is used for IXl/(N + 1) > 0.6. For IX(/(N + 1) $0.6 equation (7) is replaced by the following expansion:
I = LLIM GNW Xl = exp(- Xl c X’ tN +y!+ lj! (14)
8-n
where the value of LLIM fulfils the condition:
2 X’ N! cc ,=LLrM+ 1 W+t + 111 c,=~,c,+~
c EPSGN. (15)
This expansion allows to avoid serious numerical errors in the case when 1x1) 1 but 1X1/N is small. Equation (7) cannot bc used for IXl/(N + 1) d 0.6 because it leads to the calculation of a small difference between two large quantities.
Function HN
Function calculates HN(NS, X) from equation (6).
DISCUSSION
The program DEMO playing the role of main program, the test program TESTSL and two test subroutines, ROOT and SUFF are additional seg- ments supplied. The subroutine ROOT, uses the results of Roothaan (1951) and the subroutine SUFF uses the results of SulTczyriski (1956). Two kinds of formulas from Roothaan (1951) and Suffczyriski (1956) are used: for equal and for different exponents. The accuracy of the latter is not the best when the values of the exponents are close to each other (this does not affect the accuracy of the calculation by the present method), therefore testing in this case is not recommended.
Program availability--The program is available on request from the first author (R. Wojnecki) and can be sent by mail or e-mail (e-mail address: wojne aPLANIF61.BITNET or [emailprotected]).
REFERENCES
Davenport J. W. (1984) Phys. Rev. B 29, 2894. Femandez Rico J., Lopez R. & Ramirez G. (1988) ht.
J. Quant. Chem. 34, 121. Guseinov I. I. (1988) Phys. Rev. A 37, 2314. Homeier H. H. H. & Steinborn E. 0. (1992) Int. J. Quan?.
Chem. 41. 399. Jones H. i. (1988) Phys. Rev. A 38, 1065. Jones H. W. (1992) ht. J. Quanr. Chem. 41, 749. Kotani M., Amemiya A., Ishiguro E. & Kimura T. (1955)
Table of Molecular Integrals. Maruzen, Tokyo. Modrak P. & Wojneck; R. (1987) Phys. Rev. B 36,
5830. Modrak P. & Wojnecki R. (1989) Phys. Stat. Sol. 152,
203. Roothaan C. C. J. (1951) J. Chem. Phys. 19, 1445. Roothaan C. C. J. 0956) J. Chem. Phvs. 24, 947. Rude&erg IL (195i) J. khem. Phys. i9, 1459. Rude&erg K., Roothaan C. C. J. & Jaunzenis W. (1956)
J. Chem. Phys. 24, 201. Sharma R. R. (1976). Phys. Rev. A 13, 517. Suffczvfiski M. (1956) Acta Phvs. Pd. 15, 287. Wats& R., Da;enpdrt J. W. & Weinert M. (1986) Phys.
Rev. B 34, 8421. Watson R.. DavenDOrt
Rev. B 35, 9284.‘ J. W. & Weinert M. (1987) Phys.
Weinert M., Davenport J. W. & Watson R. E. (1986) Phys. Rev. B 34, 2971.
