1      Subroutine hf3mkr(Axyz,Bxyz,Cxyz,alpha,Gxyz,
2     &    RS,GC,ff,R,R0,IJK,Nabc,Lg,Lg3)
3c
4c $Id$
5c
6      Implicit none
7c::passed
8      integer Nabc, Lg, Lg3
9c--> Cartesian Coordinates for centers a, b, c
10
11      Double Precision Axyz(3),Bxyz(3),Cxyz(3) ! [input]
12
13c--> Exponents (1:3,*) and ES prefactors (4:4,*)
14
15      Double Precision alpha(4,Nabc) ! [input]
16
17c--> Auxiliary Function Integrals & Index
18
19      Double Precision R0(Nabc,Lg3) ! [output]
20      Integer IJK(0:Lg,0:Lg,0:Lg) ! [output]
21
22c--> Scratch Space
23
24      Double Precision Gxyz(3,Nabc), GC(Nabc,3)
25      Double Precision RS(Nabc), ff(2,Nabc), R(Nabc,0:Lg,Lg3)
26c::local
27      double precision PI, PI4
28      Parameter (PI=3.1415926535898D0,PI4=4.D0/PI)
29c
30      double precision a, b, c, abci
31*rak: double precision ab, abi
32      double precision Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz
33*rak:      Double Precision Px, Py, Pz
34      double precision GCx, GCy, GCz
35      double precision alpha_t
36*acc_debug:      double precision accy
37*acc_debug:      integer accy_cnt
38*acc_debug:      logical reached
39      integer mp, j, m, n
40
41c
42c Define the auxiliary function integrals necessary to compute the three
43c center nuclear attraction integrals (NAIs). These integrals are scaled
44c by an appropriate factor, RS, defined as
45c
46c         / a + b + c \ 1/2
47c   RS = | ----------- |
48c         \    PI/4   /
49c
50c
51c******************************************************************************
52*      call dfill((Nabc*(Lg+1)*Lg3), 0.0d00 , R, 1)
53*      call dfill((Nabc*Lg3),0.0d00, r0, 1)
54
55c Define the center "P" plus C to get "G" center.
56
57      Ax = Axyz(1)
58      Ay = Axyz(2)
59      Az = Axyz(3)
60
61      Bx = Bxyz(1)
62      By = Bxyz(2)
63      Bz = Bxyz(3)
64
65      Cx = Cxyz(1)
66      Cy = Cxyz(2)
67      Cz = Cxyz(3)
68
69      do 00100 mp = 1,Nabc
70
71        a = alpha(1,mp)
72        b = alpha(2,mp)
73        c = alpha(3,mp)
74
75*rak:        ab = a + b
76*rak:        abi = 1/ab
77*rak:
78*rak:        px = abi*(a*Ax + b*Bx)
79*rak:        py = abi*(a*Ay + b*By)
80*rak:        pz = abi*(a*Az + b*Bz)
81*rak:   abci = 1/(ab+c)
82
83        abci = 1/(a+b+c)
84
85*rak:        Gxyz(1,mp) = abci*(ab*px + c*Cx)
86*rak:        Gxyz(2,mp) = abci*(ab*py + c*Cy)
87*rak:        Gxyz(3,mp) = abci*(ab*pz + c*Cz)
88
89        Gxyz(1,mp) = abci*(a*Ax + b*Bx + c*Cx)
90        Gxyz(2,mp) = abci*(a*Ay + b*By + c*Cy)
91        Gxyz(3,mp) = abci*(a*Az + b*Bz + c*Cz)
92
93c Define the scaling factor.
94
95        RS(mp) = sqrt((a+b+c)*PI4)
96
9700100 continue
98
99c Define factors necessary to compute incomplete gamma function and the
100c auxiliary functions.
101
102       do 00200 m = 1,Nabc
103
104         alpha_t = alpha(1,m) + alpha(2,m) + alpha(3,m)
105
106         ff(1,m) = RS(m)
107         ff(2,m) = -2.D0*alpha_t
108
109         GCx = Gxyz(1,m) - Cx
110         GCy = Gxyz(2,m) - Cy
111         GCz = Gxyz(3,m) - Cz
112
113         R(m,0,1) = alpha_t*(GCx*GCx + GCy*GCy + GCz*GCz)
114
115         GC(m,1) = GCx
116         GC(m,2) = GCy
117         GC(m,3) = GCz
118
11900200  continue
120
121c Evaluate the incomplete gamma function.
122
123       call igamma(R,Nabc,Lg)
124
125*acc_debug:       accy = 1.0d-30
126*acc_debug:       accy_cnt = 0
127*acc_debug:       reached = .false.
128*acc_debug:00001  continue
129*acc_debug:       call igamma_acc(R,Nabc,Lg,accy,reached)
130*acc_debug:       if (.not.reached) then
131*acc_debug:         accy_cnt = accy_cnt + 1
132*acc_debug:         accy = accy/5.0d00
133*acc_debug:         write(6,*)' accy reset to ',accy,' count is ',accy_cnt
134*acc_debug:         goto 00001
135*acc_debug:       endif
136
137c Define the initial auxiliary functions (i.e., R000j, j=1,Lg).
138
139       do 00300 j = 0,Lg
140         do 00400 m = 1,Nabc
141           R(m,j,1) = ff(1,m)*R(m,j,1)
142           ff(1,m) = ff(1,m)*ff(2,m)
14300400    continue
14400300  continue
145
146c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0).
147
148       call hfmkr(R,IJK,GC,Nabc,Lg,Lg3)
149
150c Transfer to R0 array.
151
152       do 00500 n = 1,Lg3
153         do 00600 m = 1,Nabc
154           R0(m,n) = R(m,0,n)
15500600    continue
15600500  continue
157
158       end
159