1      Subroutine hfnai(E,R0,IJK,Vab,Nints,NPP,La,Lb,Li,Lp,Lp3,canAB)
2c $Id$
3
4      Implicit none
5
6      integer Nints,NPP,La,Lb,Li,Lp,Lp3
7      Logical canAB
8
9c--> Hermite Linear Expansion Coefficients
10
11      double precision E(3,NPP,0:((La+Li)+(Lb+Li)),0:(La+Li),0:(Lb+Li))
12
13c--> Auxiliary Function Integrals & Index
14
15      double precision R0(NPP,Lp3)
16      integer IJK(0:Lp,0:Lp,0:Lp)
17
18c--> Nuclear Attraction Integrals
19
20      double precision Vab(Nints)
21
22c--> Scratch Space
23
24      integer Nxyz(3)
25
26c--> Local variables
27
28      integer nn,ma,mb,mb_limit,mp,np,La2,Lb2
29      integer Ia,Ja,Ka, Ib,Jb,Kb, Ip,Jp,Kp
30      double precision vab_int
31c
32c Compute the nuclear attraction integrals.
33c
34c     formula:
35c           __
36c           \    Ia,Ib    Ja,Jb    Ka,Kb
37c     Vab = /  Ex     * Ey     * Ez     * R
38c           --   Ip       Jp       Kp      Ip,Jp,Kp
39c        Ip,Jp,Kp
40c
41c******************************************************************************
42
43c Initialize the block of NAIs.
44
45!      call dcopy(Nints,0d0,0,Vab,1)
46
47c Define the number of shell components on each center.
48
49      La2 = ((La+1)*(La+2))/2
50      Lb2 = ((Lb+1)*(Lb+2))/2
51
52c Loop over shell components.
53
54      nn = 0
55
56      do ma = 1,La2
57
58c Define the angular momentum indices for shell "A".
59
60        call getNxyz(La,ma,Nxyz)
61
62        Ia = Nxyz(1)
63        Ja = Nxyz(2)
64        Ka = Nxyz(3)
65
66        if( canAB )then
67          mb_limit = ma
68        else
69          mb_limit = Lb2
70        end if
71
72        do mb = 1,mb_limit
73
74c Define the angular momentum indices for shell "B".
75
76          call getNxyz(Lb,mb,Nxyz)
77
78          Ib = Nxyz(1)
79          Jb = Nxyz(2)
80          Kb = Nxyz(3)
81
82          nn = nn + 1
83
84          vab_int = 0.0d00
85
86          do Ip = 0,Ia+Ib
87            do Jp = 0,Ja+Jb
88              do Kp = 0,Ka+Kb
89
90                np = IJK(Ip,Jp,Kp)
91
92                do mp = 1,NPP
93                  vab_int = vab_int +
94     &                E(1,mp,Ip,Ia,Ib)*
95     &                E(2,mp,Jp,Ja,Jb)*
96     &                E(3,mp,Kp,Ka,Kb)*R0(mp,np)
97
98                end do
99
100              end do
101            end do
102          end do
103
104          Vab(nn) = vab_int
105
106        end do
107
108      end do
109
110      end
111*----------------------------------------------------------------------
112      Subroutine hfnai_gc(E,R0,IJK,Vab,VabP,VabH,Acoefs,Bcoefs,ipairp,
113     &    NPA,NPB,NCA,NCB,NPP,La,Lb,La2,Lb2,Li,Lp,Lp3,canAB)
114
115      implicit none
116
117      logical canAB
118      integer NPA,NPB,NCA,NCB,NPP,La,Lb,La2,Lb2,Li,Lp,Lp3
119
120c--> Hermite Linear Expansion Coefficients
121
122      double precision E(3,NPP,0:((La+Li)+(Lb+Li)),0:(La+Li),0:(Lb+Li))
123
124c--> Index of primitives
125
126      integer ipairp(2,NPP)
127
128c--> Auxiliary Function Integrals & Index
129
130      double precision R0(NPP,Lp3)
131      integer IJK(0:Lp,0:Lp,0:Lp)
132
133c--> Nuclear Attraction Integrals
134
135      double precision Vab(Lb2,NCB,La2,NCA)
136      double precision VabP(NPP)
137      double precision VabH(NPA,NCB)
138
139c--> general contraction matrices
140
141      double precision Acoefs(NPA,NCA)
142      double precision Bcoefs(NPB,NCB)
143
144c--> Scratch Space
145
146      integer Nxyz(3)
147
148c--> Local variables
149
150      integer ma,mb,mp,np, ica,icb,icb_limit,ipa,ipb
151      integer Ia,Ja,Ka, Ib,Jb,Kb, Ip,Jp,Kp
152c
153c Compute the nuclear attraction integrals.
154c
155c     formula:
156c           __
157c           \    Ia,Ib    Ja,Jb    Ka,Kb
158c     Vab = /  Ex     * Ey     * Ez     * R
159c           --   Ip       Jp       Kp      Ip,Jp,Kp
160c        Ip,Jp,Kp
161c
162c******************************************************************************
163
164c Initialize the block of NAIs.
165
166      call dcopy(La2*Lb2*nca*ncb,0d0,0,Vab,1)
167
168c Loop over shell components.
169
170      do ma = 1,La2
171
172c Define the angular momentum indices for shell "A".
173
174        call getNxyz(La,ma,Nxyz)
175
176        Ia = Nxyz(1)
177        Ja = Nxyz(2)
178        Ka = Nxyz(3)
179
180        do mb = 1,Lb2
181
182c Define the angular momentum indices for shell "B".
183
184          call getNxyz(Lb,mb,Nxyz)
185
186          Ib = Nxyz(1)
187          Jb = Nxyz(2)
188          Kb = Nxyz(3)
189          call dcopy(NPP,0d0,0,VabP,1)
190          do Ip = 0,Ia+Ib
191            do Jp = 0,Ja+Jb
192              do Kp = 0,Ka+Kb
193
194                np = IJK(Ip,Jp,Kp)
195
196                do mp = 1,NPP
197                  VabP(mp) = VabP(mp) + R0(mp,np)*
198     &                E(1,mp,Ip,Ia,Ib)*
199     &                E(2,mp,Jp,Ja,Jb)*
200     &                E(3,mp,Kp,Ka,Kb)
201                end do
202
203              end do
204            end do
205          end do
206
207c Contract over B shell
208
209          call dcopy(NCB*NPA,0d0,0,VabH,1)
210          do icb = 1,NCB
211            do mp = 1,NPP
212              ipa = ipairp(1,mp)
213              ipb = ipairp(2,mp)
214              VabH(ipa,icb) = VabH(ipa,icb)+VabP(mp)*Bcoefs(ipb,icb)
215            end do
216          end do
217
218c Contract over A shell
219
220          do ica = 1,NCA
221            if( canAB )then
222              icb_limit = ica
223            else
224              icb_limit = NCB
225            end if
226            do icb = 1,icb_limit
227              do ipa = 1,NPA
228                Vab(mb,icb,ma,ica) = Vab(mb,icb,ma,ica) +
229     &              VabH(ipa,icb)*Acoefs(ipa,ica)
230              end do
231            end do
232          end do
233
234        end do
235
236      end do
237
238      if (canAB) call canon_ab(Vab,Vab,Lb2*NCB,La2*NCA)
239
240      end
241************************************************************************
242      subroutine canon_ab (Vcanon,Vsquare,NB,NA)
243      implicit none
244      integer NA,NB,i,j,k
245      double precision Vcanon(*),Vsquare(NB,NA)
246
247      k = 1
248      do i = 1,NA
249        do j = 1,i
250          k = k+1
251          Vcanon(k) = Vsquare(j,i)
252        end do
253      end do
254
255      end
256