1C$PRAGMA SUN OPT=2 2 Subroutine hf1mke3(Axyz,Bxyz,Cxyz,alpha,G,GT,ABC2I,E, 3 & NABC,La,Lb,Lc) 4c $Id$ 5 6 Implicit real*8 (a-h,o-z) 7 Implicit integer (i-n) 8 9 Dimension Axyz(3),Bxyz(3),Cxyz(3),alpha(4,NABC),G(3,NABC) 10 11 Dimension GT(NABC,3),ABC2I(3*NABC) 12 13 Dimension E(NABC,3,0:(La+Lb+Lc),0:La,0:Lb,0:Lc) 14c 15c Define the Hermite linear expansion coefficients for the product of 16c three monomials. These coefficients are scaled by a factor appropriate 17c for 3-ctr OIs, defined as 18c 19c / PI \ 3/2 / a b __2 \ / (a + b) c __2 \ 20c ES = | ----------- | EXP| - ----- AB | EXP| - ----------- PC | 21c \ a + b + c / \ a + b / \ a + b + c / 22c 23c This scaling factor is passed as the 4th index of the exponents array 24c (i.e., "alpha(4,m)"). 25c 26c Recursion Formulae: 27c 28c 0,0,0 29c Ex = ES 30c 0 31c 32c 0,0,0 33c Ey = 1.0 34c 0 35c 36c 0,0,0 37c Ez = 1.0 38c 0 39c 40c Ia+1,Ib,Ic Ia,Ib,Ic Ia,Ib,Ic Ia,Ib,Ic 41c Ex = ABC2I Ex + GAx Ex + (Ip+1) Ex 42c Ip Ip-1 Ip Ip+1 43c 44c Ia,Ib+1,Ic Ia,Ib,Ic Ia,Ib,Ic Ia,Ib,Ic 45c Ex = ABC2I Ex + GBx Ex + (Ip+1) Ex 46c Ip Ip-1 Ip Ip+1 47c 48c Ia,Ib,Ic+1 Ia,Ib,Ic Ia,Ib,Ic Ia,Ib,Ic 49c Ex = ABC2I Ex + GCx Ex + (Ip+1) Ex 50c Ip Ip-1 Ip Ip+1 51c 52c Indices for E(m,k,Ip,Ia,Ib,Ic): 53c 54c 1 [ m ] NABC 55c 1 [ k ] 3 {X,Y,Z} 56c 0 [ Ip ] Ia+Ib+Ic 57c 0 [ Ia ] La 58c 0 [ Ib ] Lb 59c 0 [ Ic ] Lc 60c 61c****************************************************************************** 62 63c Initialize the Hermite expansion coefficients. 64 65 isz_e = NABC*3*(La+Lb+Lc+1)*(La+1)*(Lb+1)*(Lc+1) 66 call dcopy(isz_e,0d0,0,E,1) 67c Define E(Ip,Ia,Ib,Ic) for Ip=0, Ia=0, Ib=0, Ic=0. 68 69 do 100 m = 1,NABC 70 71 E(m,1,0,0,0,0) = alpha(4,m) 72 E(m,2,0,0,0,0) = 1.D0 73 E(m,3,0,0,0,0) = 1.D0 74 75 ABC2I(m ) = 0.5D0/(alpha(1,m) + alpha(2,m) + alpha(3,m)) 76 ABC2I(m+1*NABC) = ABC2I(m) 77 ABC2I(m+2*NABC) = ABC2I(m) 78 79 GT(m,1) = G(1,m) - Axyz(1) 80 GT(m,2) = G(2,m) - Axyz(2) 81 GT(m,3) = G(3,m) - Axyz(3) 82 83 100 continue 84 85c Define E(Ip,Ia,Ib,Ic) for Ip=0,Ia+Ib+Ic, Ia=1,La, Ib=0, Ic=0. 86 87 do 250 Ia = 1,La 88 89c ===> Ip = 0 <=== 90 91 do 210 m = 1,3*NABC 92 E(m,1,0,Ia,0,0) = GT(m,1)*E(m,1,0,Ia-1,0,0) 93 & + E(m,1,1,Ia-1,0,0) 94 210 continue 95 96c ===> Ip = 1,Ia-1 <=== 97 98 do 230 Ip = 1,Ia-1 99 100 do 220 m = 1,3*NABC 101 E(m,1,Ip,Ia,0,0) = ABC2I(m)*E(m,1,Ip-1,Ia-1,0,0) 102 & + GT(m,1)*E(m,1,Ip ,Ia-1,0,0) 103 & + (Ip+1)*E(m,1,Ip+1,Ia-1,0,0) 104 220 continue 105 106 230 continue 107 108c ===> Ip = Ia <=== 109 110 Ip = Ia 111 112 do 240 m = 1,3*NABC 113 E(m,1,Ip,Ia,0,0) = ABC2I(m)*E(m,1,Ip-1,Ia-1,0,0) 114 & + GT(m,1)*E(m,1,Ip ,Ia-1,0,0) 115 240 continue 116 117 250 continue 118 119c Define E(Ip,Ia,Ib,Ic) for Ip=0,Ia+Ib+Ic, Ia=0,La, Ib=1,Lb, Ic=0. 120 121 do 300 m = 1,NABC 122 GT(m,1) = G(1,m) - Bxyz(1) 123 GT(m,2) = G(2,m) - Bxyz(2) 124 GT(m,3) = G(3,m) - Bxyz(3) 125 300 continue 126 127 do 360 Ib = 1,Lb 128 129 do 350 Ia = 0,La 130 131c ===> Ip = 0 <=== 132 133 do 310 m = 1,3*NABC 134 E(m,1,0,Ia,Ib,0) = GT(m,1)*E(m,1,0,Ia,Ib-1,0) 135 & + E(m,1,1,Ia,Ib-1,0) 136 310 continue 137 138c ===> Ip = 1,Ia+Ib-1 <=== 139 140 do 330 Ip = 1,Ia+Ib-1 141 142 do 320 m = 1,3*NABC 143 E(m,1,Ip,Ia,Ib,0) = ABC2I(m)*E(m,1,Ip-1,Ia,Ib-1,0) 144 & + GT(m,1)*E(m,1,Ip ,Ia,Ib-1,0) 145 & + (Ip+1)*E(m,1,Ip+1,Ia,Ib-1,0) 146 320 continue 147 148 330 continue 149 150c ===> Ip = Ia+Ib <=== 151 152 Ip = Ia + Ib 153 154 do 340 m = 1,3*NABC 155 E(m,1,Ip,Ia,Ib,0) = ABC2I(m)*E(m,1,Ip-1,Ia,Ib-1,0) 156 & + GT(m,1)*E(m,1,Ip ,Ia,Ib-1,0) 157 340 continue 158 159 350 continue 160 161 360 continue 162 163c Define E(Ip,Ia,Ib,Ic) for Ip=0,Ia+Ib+Ic, Ia=0,La, Ib=0,Lb, Ic=1,Lc. 164 165 do 400 m = 1,NABC 166 GT(m,1) = G(1,m) - Cxyz(1) 167 GT(m,2) = G(2,m) - Cxyz(2) 168 GT(m,3) = G(3,m) - Cxyz(3) 169 400 continue 170 171 do 470 Ic = 1,Lc 172 173 do 460 Ib = 0,Lb 174 175 do 450 Ia = 0,La 176 177c ===> Ip = 0 <=== 178 179 do 410 m = 1,3*NABC 180 E(m,1,0,Ia,Ib,Ic) = GT(m,1)*E(m,1,0,Ia,Ib,Ic-1) 181 & + E(m,1,1,Ia,Ib,Ic-1) 182 410 continue 183 184c ===> Ip = 1,Ia+Ib+Ic-1 <=== 185 186 do 430 Ip = 1,Ia+Ib+Ic-1 187 188 do 420 m = 1,3*NABC 189 E(m,1,Ip,Ia,Ib,Ic) = ABC2I(m)*E(m,1,Ip-1,Ia,Ib,Ic-1) 190 & + GT(m,1)*E(m,1,Ip ,Ia,Ib,Ic-1) 191 & + (Ip+1)*E(m,1,Ip+1,Ia,Ib,Ic-1) 192 420 continue 193 194 430 continue 195 196c ===> Ip = Ia+Ib+Ic <=== 197 198 Ip = Ia + Ib + Ic 199 200 do 440 m = 1,3*NABC 201 E(m,1,Ip,Ia,Ib,Ic) = ABC2I(m)*E(m,1,Ip-1,Ia,Ib,Ic-1) 202 & + GT(m,1)*E(m,1,Ip ,Ia,Ib,Ic-1) 203 440 continue 204 205 450 continue 206 207 460 continue 208 209 470 continue 210 211 end 212