1      Subroutine hfset(Axyz,Aprims,Acoefs,NPA,NCA,
2     &                 Bxyz,Bprims,Bcoefs,NPB,NCB,
3     &                 GENCON,alpha,ipair,ES,NPP)
4! $Id$
5
6      Implicit real*8 (a-h,o-z)
7      Implicit integer (i-n)
8
9      Logical GENCON
10
11      Parameter (PI=3.1415926535898d0)
12      Parameter (EXPLIM=100.d0)
13#include "apiP.fh"
14
15!--> Cartesian Coordinates, Primitives & Contraction Coefficients
16
17      Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
18      Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
19
20!--> Exponents, Pair Index & Prefactors for 2-ctr Overlap Distributions
21
22      Dimension alpha(2,(NPA*NPB)),ipair(2,(NPA*NPB)),ES(3,(NPA*NPB))
23!
24! Compute the prefactors of the overlap distributions formed by the product of
25! two primitive Gaussians. These prefactors are defined as
26!
27!   ES = ESx * ESy * ESz, where
28!
29!
30!          /  PI   \ 1/2      /   a b     2 \
31!   ESx = | ------- |     EXP| - -----  Rx   |
32!          \ a + b /          \  a + b      /
33!
34!
35! N.B. 1) Overlap distributions with prefactors less than a given tolerance
36!         are removed from the list. This shortened list is of length "NPP".
37!      2) For segmented contractions, the product of contraction coefficients
38!         is also incorporated in the prefactor.
39!
40!******************************************************************************
41
42      ! Python: 1.5*math.log(math.pi) = 1.7170948287741004
43      parameter (LOGPIx15=1.7170948287741d0)
44      parameter (N3OVER2=-1.5d0)
45      double precision eps_small
46      parameter (eps_small=1d-32)
47      double precision const
48      if(val_int_acc.lt.eps_small)  then
49         const = LOGPIx15 - log(eps_small)
50      else
51         const = LOGPIx15 - log(val_int_acc)
52      endif
53
54
55      Rx2 = (Axyz(1) - Bxyz(1))**2
56      Ry2 = (Axyz(2) - Bxyz(2))**2
57      Rz2 = (Axyz(3) - Bxyz(3))**2
58
59#if 0
60
61      ! ORIGINAL CODE BELOW
62
63      m2 = 0
64      do mpa = 1,NPA
65        a = Aprims(mpa)
66        do mpb = 1,NPB
67          b = Bprims(mpb)
68
69          abi = 1.0d0/(a+b)
70          beta = a*b*abi
71
72          if ( ( beta*(Rx2+Ry2+Rz2) ) .lt.
73     &         ( const + N3OVER2 * log(a+b) ) ) then
74
75            m2 = m2 + 1
76
77            alpha(1,m2) = a
78            alpha(2,m2) = b
79
80            ipair(1,m2) = mpa
81            ipair(2,m2) = mpb
82
83            s = sqrt(PI*abi)
84
85            ESx = s*exp(-min(EXPLIM,beta*Rx2))
86            ESy = s*exp(-min(EXPLIM,beta*Ry2))
87            ESz = s*exp(-min(EXPLIM,beta*Rz2))
88
89            ! GENCON is an easy branch to predict but should
90            ! still be hoisted out of the loops.
91            if( GENCON )then
92              ES(1,m2) = ESx
93              ES(2,m2) = ESy
94              ES(3,m2) = ESz
95            else ! GENCON
96              ES(1,m2) = ESx*(Acoefs(mpa,1)*Bcoefs(mpb,1))
97              ES(2,m2) = ESy
98              ES(3,m2) = ESz
99            end if ! GENCON
100
101          end if ! threshold
102
103        enddo ! mpb
104      enddo ! mpa
105      NPP = m2
106
107#else
108
109      m2 = 0
110      do mpa = 1,NPA
111        a = Aprims(mpa)
112        do mpb = 1,NPB
113          b = Bprims(mpb)
114
115          abi = 1.0d0/(a+b)
116          beta = a*b*abi
117          if ( ( beta*(Rx2+Ry2+Rz2) ) .lt.
118     &         ( const + N3OVER2 * log(a+b) ) ) then
119
120            m2 = m2 + 1
121
122            alpha(1,m2) = a
123            alpha(2,m2) = b
124
125            ipair(1,m2) = mpa
126            ipair(2,m2) = mpb
127
128          end if ! threshold
129        enddo ! mpb
130      enddo ! mpa
131      NPP = m2
132
133      if (GENCON) then
134        do m2 = 1, NPP
135
136          a = alpha(1,m2)
137          b = alpha(2,m2)
138
139          abi = 1.0d0/(a+b)
140          beta = a*b*abi
141
142          s = sqrt(PI*abi)
143
144          ES(1,m2) = s*exp(-min(EXPLIM,beta*Rx2))
145          ES(2,m2) = s*exp(-min(EXPLIM,beta*Ry2))
146          ES(3,m2) = s*exp(-min(EXPLIM,beta*Rz2))
147
148        enddo ! m2
149      else ! GENCON
150        do m2 = 1, NPP
151
152          a = alpha(1,m2)
153          b = alpha(2,m2)
154
155          abi = 1.0d0/(a+b)
156          beta = a*b*abi
157
158          mpa = ipair(1,m2)
159          mpb = ipair(2,m2)
160
161          s = sqrt(PI*abi)
162
163          c = Acoefs(mpa,1)
164          d = Bcoefs(mpb,1)
165
166          ESx = s*exp(-min(EXPLIM,beta*Rx2))
167          ESy = s*exp(-min(EXPLIM,beta*Ry2))
168          ESz = s*exp(-min(EXPLIM,beta*Rz2))
169
170          !ES(1,m2) = ESx*(Acoefs(mpa,1)*Bcoefs(mpb,1))
171          ES(1,m2) = ESx*c*d
172          ES(2,m2) = ESy
173          ES(3,m2) = ESz
174
175        enddo ! m2
176      endif ! GENCON
177
178#endif
179
180!      write(6,*)'-----start------ pair matrix '
181!      do i=1,NPP
182!        write(6,*)i,ipair(1,i),ipair(2,i)
183!      enddo
184!      write(6,*)'----- end ------ pair matrix '
185      end
186