1 Subroutine hf1set3(Axyz,Aprims,Acoef,NPA, 2 & Bxyz,Bprims,Bcoef,NPB, 3 & Cxyz,Cprims,Ccoef,NPC, 4 & alpha,NABC) 5c $Id$ 6 7 Implicit real*8 (a-h,o-z) 8 Implicit integer (i-n) 9 10 Parameter (PI=3.1415926535898D0) 11 Parameter (EXPLIM=100.D0) 12#include "apiP.fh" 13 14 Dimension Axyz(3),Aprims(NPA),Acoef(NPA) 15 Dimension Bxyz(3),Bprims(NPB),Bcoef(NPB) 16 Dimension Cxyz(3),Cprims(NPC),Ccoef(NPC) 17 18 Dimension alpha(4,NPA*NPB*NPC) 19c 20c The prefactor of the charge distribution, formed by the product of 21c three Gaussians, is defined as 22c 23c / PI \ 3/2 / a b __2 \ / (a + b) c __2 \ 24c ES = | ----------- | EXP| - ----- AB | EXP| - ----------- PC | 25c \ a + b + c / \ a + b / \ a + b + c / 26c 27c N.B. 1) This prefactor is returned as the 4th index of the exponents array. 28c (i.e., "alpha(4,m)") 29c 2) Charge distributions, whose prefactor is less than a given tolerance, 30c are removed from the list. The shortened list is of length "NABC". 31c 3) The product of contraction coefficients is also incorporated in 32c the prefactor. For a general contraction implementation, this 33c should not be done at this point. 34c 35c****************************************************************************** 36 37 m = 0 38 do 10 ipa = 1,NPA 39 do 10 ipb = 1,NPB 40 do 10 ipc = 1,NPC 41 m = m + 1 42 43 alpha(1,m) = Aprims(ipa) 44 alpha(2,m) = Bprims(ipb) 45 alpha(3,m) = Cprims(ipc) 46 47 alpha(4,m) = Acoef(ipa)*Bcoef(ipb)*Ccoef(ipc) 48 49 10 continue 50 51 Ax = Axyz(1) 52 Ay = Axyz(2) 53 Az = Axyz(3) 54 55 Bx = Bxyz(1) 56 By = Bxyz(2) 57 Bz = Bxyz(3) 58 59 Cx = Cxyz(1) 60 Cy = Cxyz(2) 61 Cz = Cxyz(3) 62 63 RABsq = (Ax-Bx)**2 + (Ay-By)**2 + (Az-Bz)**2 64 65 m2 = 0 66 do 20 m1 = 1,NPA*NPB*NPC 67 68 A = alpha(1,m1) 69 B = alpha(2,m1) 70 C = alpha(3,m1) 71 72 AB = A + B 73 ABI = 1/AB 74 75 ES = exp(-min(EXPLIM,A*B*ABI*RABsq)) 76 77 Px = ABI*(A*Ax + B*Bx) 78 Py = ABI*(A*Ay + B*By) 79 Pz = ABI*(A*Az + B*Bz) 80 81 ABCI = 1/(A + B + C) 82 83 RPCsq = (Px-Cx)**2 + (Py-Cy)**2 + (Pz-Cz)**2 84 85 ES = ES*exp(-min(EXPLIM,(AB*C*ABCI)*RPCsq)) 86 87 s = sqrt(PI*ABCI) 88 89 ES = s*s*s*ES 90 91 if( ES .gt. val_int_acc )then 92 93 m2 = m2 + 1 94 95 alpha(1,m2) = A 96 alpha(2,m2) = B 97 alpha(3,m2) = C 98 99 alpha(4,m2) = ES*alpha(4,m1) 100 101 end if 102 103 20 continue 104 NABC = m2 105 106 end 107