1      Subroutine hf2mkr(Axyz,Bxyz,Cxyz,Dxyz,alpha_p,alpha_q,R0,IJK,
2     &    Pxyz,Qxyz,PQ,ff,R,ffscr,rscr,NPP,NPQ,Lr,Lr3)
3c
4c $Id$
5
6      implicit none
7#include "case.fh"
8c
9      double precision PI,P1
10      Parameter (PI=3.1415926535898D0,P1=4.D0/PI)
11
12c--> ints
13      integer npp,npq
14      integer Lr, Lr3
15      integer IJK(0:Lr,0:Lr,0:Lr)
16
17c--> Cartesian Coordinates
18      double precision Axyz(3),Bxyz(3),Cxyz(3),Dxyz(3)
19
20c--> Exponents
21
22      double precision alpha_p(2,NPP),alpha_q(2,NPQ)
23
24c--> Auxiliary Function Integrals & Index
25
26      double precision R0((NPP*NPQ),Lr3)
27
28c--> Scratch Space
29
30      double precision Pxyz(3,NPP),Qxyz(3,NPQ),PQ((NPP*NPQ),3),
31     &  ff(2,(NPP*NPQ)), R((NPP*NPQ),0:Lr,Lr3)
32c
33      double precision ffscr(2,(NPP*NPQ)), Rscr((NPP*NPQ),0:Lr,Lr3),
34     & rhoscaled,radd
35
36      integer mp, mq, mr, j, n
37      double precision a,b,c,d,f1,f2,p,q,rho
38      double precision PQx, PQy, PQz
39c
40c Evaluate the auxiliary function integrals. These integrals are scaled by a
41c factor appropriate for ERIs, defined as
42c
43c                          1/2
44c  ES = ((4/PI)*(pq/(p+q)))           where p = a + b and q = c + d.
45c
46c******************************************************************************
47c
48c Define the center "P".
49
50      do 100 mp = 1,NPP
51
52       a = alpha_p(1,mp)
53       b = alpha_p(2,mp)
54
55       f1 = a/(a+b)
56       f2 = b/(a+b)
57
58       Pxyz(1,mp) = f1*Axyz(1) + f2*Bxyz(1)
59       Pxyz(2,mp) = f1*Axyz(2) + f2*Bxyz(2)
60       Pxyz(3,mp) = f1*Axyz(3) + f2*Bxyz(3)
61
62  100 continue
63
64c Define the center "Q".
65
66      do 110 mq = 1,NPQ
67
68       c = alpha_q(1,mq)
69       d = alpha_q(2,mq)
70
71       f1 = c/(c+d)
72       f2 = d/(c+d)
73
74       Qxyz(1,mq) = f1*Cxyz(1) + f2*Dxyz(1)
75       Qxyz(2,mq) = f1*Cxyz(2) + f2*Dxyz(2)
76       Qxyz(3,mq) = f1*Cxyz(3) + f2*Dxyz(3)
77
78  110 continue
79
80c Define factors necessary to compute incomplete gamma function and the
81c auxiliary functions.
82
83      mr = 0
84      do 125 mp = 1,NPP
85        do 120 mq = 1,NPQ
86          mr = mr + 1
87
88          p = alpha_p(1,mp) + alpha_p(2,mp)
89          q = alpha_q(1,mq) + alpha_q(2,mq)
90
91          rho = p*q/(p+q)
92
93          PQx = Pxyz(1,mp) - Qxyz(1,mq)
94          PQy = Pxyz(2,mp) - Qxyz(2,mq)
95          PQz = Pxyz(3,mp) - Qxyz(3,mq)
96
97          PQ(mr,1) = PQx
98          PQ(mr,2) = PQy
99          PQ(mr,3) = PQz
100c
101          R(mr,0,1) = rho*(PQx**2 + PQy**2 + PQz**2)
102          ff(1,mr) = 2.d0*dsqrt(rho/pi)
103          ff(2,mr) = -2.D0*rho
104c
105          if (doscreen) then
106           rhoscaled = rho
107           call case_md(rhoscaled)
108           Rscr(mr,0,1) = rhoscaled*(PQx**2 + PQy**2 + PQz**2)
109           ffscr(1,mr) = 2.d0*dsqrt(rhoscaled/pi)
110           ffscr(2,mr) = -2.D0*rhoscaled
111          end if
112c
113  120   continue
114  125 continue
115
116c Evaluate the incomplete gamma function.
117
118      call igamma(R,(NPP*NPQ),Lr)
119      if (doscreen) call igamma(Rscr,(NPP*NPQ),Lr)
120
121c Define the initial auxiliary functions (i.e., R000j, j=1,Lr).
122
123      do 135 j = 0,Lr
124        do 130 mr = 1,(NPP*NPQ)
125             R(mr,j,1) = ff(1,mr)*R(mr,j,1)
126             ff(1,mr) = ff(1,mr)*ff(2,mr)
127  130   continue
128             if (doscreen) then
129                do  mr = 1,(NPP*NPQ)
130                   Rscr(mr,j,1) = ffscr(1,mr)*Rscr(mr,j,1)
131                   ffscr(1,mr) = ffscr(1,mr)*ffscr(2,mr)
132                enddo
133             endif
134  135 continue
135
136c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0).
137
138      call hfmkr(R,IJK,PQ,(NPP*NPQ),Lr,Lr3)
139      if (doscreen) call hfmkr(Rscr,IJK,PQ,(NPP*NPQ),Lr,Lr3)
140
141c Transfer to R0 array.
142
143c       write(6,*) "cam_alpha:",cam_alpha
144c       write(6,*) "cam_beta:",cam_beta
145c       write(6,*) "doscreen:",doscreen
146c       write(6,*) "cam_srhf:",cam_srhf
147c
148      if (doscreen) then
149         do  n = 1,Lr3
150            if (cam_srhf) then  ! for pure short-range HF (1-erf(r)/r)
151               do  mr = 1,(NPP*NPQ)
152                  R0(mr,n) = R(mr,0,n)-Rscr(mr,0,n)
153               enddo
154            else
155               do  mr = 1,(NPP*NPQ)
156                  R0(mr,n) = cam_alpha*R(mr,0,n)+cam_beta*Rscr(mr,0,n)
157               enddo
158            endif
159         enddo
160      else
161         do  n = 1,Lr3
162            do  mr = 1,(NPP*NPQ)
163               R0(mr,n) = R(mr,0,n)
164            enddo
165         enddo
166      endif
167
168      end
169