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