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