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