1      Subroutine hf1mkr(Axyz,Bxyz,Cxyz,zan,exinv,ncenters,
2     &                  alpha,Pxyz,RS,PC,ff,R,
3     &                  R0,R0C,IJK,NPP,Lp,Lp3,CENTER)
4c $Id$
5
6      Implicit none
7c::passed
8      Integer ncenters                  ! number of nuclei
9      Integer NPP                       ! number of primitive products
10      Integer Lp                        ! angular momentum sum La+Lb+MXD
11      Integer Lp3                       ! number of angular momentum functions
12      Double precision Axyz(3),Bxyz(3)  ! Cartesian coordinates of centers
13      Double precision Cxyz(3,ncenters) ! Nuclear cartesian coordinates
14      Double precision zan(ncenters)    ! Nuclear charges
15      Double precision exinv(ncenters)  ! Nuclear inverse exponents
16      Double precision alpha(2,NPP)     ! basis function exponents copy
17      Double precision Pxyz(3,NPP)      ! coordinates of product center
18      Double precision RS(NPP)          ! scale factor
19      Double precision PC(NPP,3)        ! distance from product center to nucleus
20      Double precision ff(2,NPP)        ! recursion factors
21      Double precision R(NPP,0:Lp,Lp3)  ! R/R000j/RIKJj functions
22      Double precision R0(NPP,Lp3)      ! R integrals
23      Double precision R0C(ncenters,NPP,Lp3) ! R integrals per center
24      integer IJK(0:Lp,0:Lp,0:Lp)       ! index array
25      Logical CENTER
26c::local
27      integer ic,j,m,n
28      Double precision a,b,f1,f2,alpha_t,beta,PCx,PCy,PCz
29      Double precision PI, PI4          ! pi, 4/pi
30      Parameter (PI=3.1415926535898D0,PI4=4.D0/PI)
31c
32c Define the auxiliary function integrals necessary to compute the nuclear
33c attraction integrals (NAIs). These integrals are scaled by an appropriate
34c factor, RS, defined as
35c
36c         / a + b \ 1/2
37c   RS = | ------- |
38c         \  PI/4 /
39c
40c The scale factor for the Hermite expansion coefficients is assumed to be
41c
42c         /   PI  \ 3/2     /   a b   __2 \
43c   ES = | ------- |    EXP| - -----  AB   |
44c         \ a + b /         \  a + b      /
45c
46c Therefore,
47c
48c            2 PI        /   a b   __2 \
49c   ES RS = -------  EXP| - -----  AB   |
50c            a + b       \  a + b      /
51c
52c******************************************************************************
53
54      do 101 m = 1,NPP
55       do 100 j = 1,Lp3
56        R0(m,j) = 0.D0
57  100  continue
58  101 continue
59
60      do 110 m = 1,NPP
61
62c Define the center "P".
63
64       a = alpha(1,m)
65       b = alpha(2,m)
66
67       f1 = a/(a+b)
68       f2 = b/(a+b)
69
70       Pxyz(1,m) = f1*Axyz(1) + f2*Bxyz(1)
71       Pxyz(2,m) = f1*Axyz(2) + f2*Bxyz(2)
72       Pxyz(3,m) = f1*Axyz(3) + f2*Bxyz(3)
73
74c Define the scaling factor.
75
76       RS(m) = sqrt((a+b)*PI4)
77
78  110 continue
79
80c Sum over all centers.
81
82      do 150 ic = 1,ncenters
83
84       beta = exinv(ic)
85
86c Define factors necessary to compute incomplete gamma function and the
87c auxiliary functions.
88
89       do 120 m = 1,NPP
90
91        alpha_t = alpha(1,m) + alpha(2,m)
92        alpha_t = alpha_t/(1.0d0+alpha_t*beta)
93
94        ff(1,m) = RS(m)
95        ff(2,m) = -2.D0*alpha_t
96
97        PCx = Pxyz(1,m) - Cxyz(1,ic)
98        PCy = Pxyz(2,m) - Cxyz(2,ic)
99        PCz = Pxyz(3,m) - Cxyz(3,ic)
100
101        R(m,0,1) = alpha_t*(PCx**2 + PCy**2 + PCz**2)
102
103        PC(m,1) = PCx
104        PC(m,2) = PCy
105        PC(m,3) = PCz
106
107  120  continue
108
109c Evaluate the incomplete gamma function.
110
111       call igamma(R,NPP,Lp)
112
113c Scale for finite nuclear size if necessary
114
115       if (beta .ne. 0.0D0) then
116        do 125 m = 1,NPP
117         alpha_t = sqrt(1.0d0+(alpha(1,m)+alpha(2,m))*beta)
118         do 126 j=0,Lp
119          R(m,j,1) = R(m,j,1)/alpha_t
120  126    continue
121  125   continue
122       end if
123
124c Define the initial auxiliary functions (i.e., R000j, j=1,Lr).
125
126       do 135 j = 0,Lp
127        do 130 m = 1,NPP
128         R(m,j,1) = ff(1,m)*R(m,j,1)
129         ff(1,m) = ff(1,m)*ff(2,m)
130  130   continue
131  135  continue
132
133c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0).
134
135       call hfmkr(R,IJK,PC,NPP,Lp,Lp3)
136
137c Transfer to R0 array.
138
139       if( CENTER )then
140        do 141 n = 1,Lp3
141         do 140 m = 1,NPP
142          R0C(ic,m,n) = -zan(ic)*R(m,0,n)
143          R0(m,n) = R0(m,n) + R0C(ic,m,n)
144  140    continue
145  141   continue
146       else
147        do 146 n = 1,Lp3
148         do 145 m = 1,NPP
149          R0(m,n) = R0(m,n) - zan(ic)*R(m,0,n)
150  145    continue
151  146   continue
152       end if
153
154  150 continue
155
156      end
157