1 Subroutine hfefi(E,R0C,IJK,bEFI, 2 & NPP,Nint,La,Lb,Li,Lp,Lp3,ncenters, 3 & MXD,canAB,ictrA,ictrB) 4c $Id$ 5 Implicit real*8 (a-h,o-z) 6 Implicit integer (i-n) 7 8 Logical canAB 9 10c--> Hermite Linear Expansion Coefficients 11 12 Dimension E(3,NPP,0:MXD,0:((La+Li)+(Lb+Li)),0:(La+Li),0:(Lb+Li)) 13 14c--> Auxiliary Function Integrals & Index 15 16 Dimension R0C(ncenters,NPP,Lp3),IJK(0:Lp,0:Lp,0:Lp) 17 18c--> Block of Electric Field Integrals 19 20 Dimension bEFI(Nint,3,ncenters) 21 22c--> Scratch Space 23 24 Dimension Nxyz(3) 25c 26c Compute the electric field integrals. 27c 28c formula: 29c __ 30c Cx \ Ia,Ib Ja,Jb Ka,Kb C 31c EFI = / Ex * Ey * Ez * R 32c ab -- Ip Jp Kp Ip+1,Jp,Kp 33c Ip,Jp,Kp 34c 35c****************************************************************************** 36 37c Initialize the block of EFIs. 38 39 do 10 ic = 1,ncenters 40 do 10 nn = 1,Nint 41 bEFI(nn,1,ic) = 0.D0 42 bEFI(nn,2,ic) = 0.D0 43 bEFI(nn,3,ic) = 0.D0 4410 continue 45 46c Define the number of shell components on each center. 47 48 La2 = ((La+1)*(La+2))/2 49 Lb2 = ((Lb+1)*(Lb+2))/2 50 51c Loop over shell components. 52 53 nn = 0 54 55 do 50 ma = 1,La2 56 57c Define the angular momentum indices for shell "A". 58 59 call getNxyz(La,ma,Nxyz) 60 61 Ia = Nxyz(1) 62 Ja = Nxyz(2) 63 Ka = Nxyz(3) 64 65 if( canAB )then 66 mb_limit = ma 67 else 68 mb_limit = Lb2 69 end if 70 71 do 40 mb = 1,mb_limit 72 73c Define the angular momentum indices for shell "B". 74 75 call getNxyz(Lb,mb,Nxyz) 76 77 Ib = Nxyz(1) 78 Jb = Nxyz(2) 79 Kb = Nxyz(3) 80 81 nn = nn + 1 82 83 do 30 Ip = 0,Ia+Ib 84 do 30 Jp = 0,Ja+Jb 85 do 30 Kp = 0,Ka+Kb 86 87 npx = IJK(Ip+1,Jp ,Kp ) 88 npy = IJK(Ip ,Jp+1,Kp ) 89 npz = IJK(Ip ,Jp ,Kp+1) 90 91 do 20 mp = 1,NPP 92 93 E3 = E(1,mp,0,Ip,Ia,Ib)* 94 & E(2,mp,0,Jp,Ja,Jb)* 95 & E(3,mp,0,Kp,Ka,Kb) 96 97 do 15 ic = 1,ncenters 98 bEFI(nn,1,ic) = bEFI(nn,1,ic) - E3*R0C(ic,mp,npx) 99 bEFI(nn,2,ic) = bEFI(nn,2,ic) - E3*R0C(ic,mp,npy) 100 bEFI(nn,3,ic) = bEFI(nn,3,ic) - E3*R0C(ic,mp,npz) 10115 continue 102c 103* bEFI(nn,1,ictra) = bEFI(nn,1,ictra) - E3*R0C(ictra,mp,npx) 104* bEFI(nn,2,ictra) = bEFI(nn,2,ictra) - E3*R0C(ictra,mp,npy) 105* bEFI(nn,3,ictra) = bEFI(nn,3,ictra) - E3*R0C(ictra,mp,npz) 106* if (ictra.ne.ictrb) then 107* bEFI(nn,1,ictrb) = bEFI(nn,1,ictrb) - E3*R0C(ictrb,mp,npx) 108* bEFI(nn,2,ictrb) = bEFI(nn,2,ictrb) - E3*R0C(ictrb,mp,npy) 109* bEFI(nn,3,ictrb) = bEFI(nn,3,ictrb) - E3*R0C(ictrb,mp,npz) 110* endif 11120 continue 112 11330 continue 114 11540 continue 116 11750 continue 118 119 end 120