1 Subroutine hf1d_cosmo( 2 & Axyz,Aprims,Acoefs,NPA,NCA,La,ictrA, 3 & Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,ictrB, 4 & Cxyz,zan, 5 & bNAI,Nint, 6 & canAB,dryrun,W0,maxW0) 7c $Id$ 8c 9c Calculate derivatives of an nuclear attraction type integral. 10c The integral is: 11c 12c E = int(dr;Chi_A(r)Chi_B(r)Q_C/|r-R_C|) 13c 14c We are interested in the derivatives with respect to the centers 15c A, B, and C. So upon return we want 16c 17c bNAI(1:nint,1:3) = dE/dA (1) 18c 19c bNAI(1:nint,4:6) = dE/dB (2) 20c 21c bNAI(1:nint,7:9) = dE/dC = -dE/dA-dE/dB (3) 22c 23c Typically the calculation of derivative integrals relies on the 24c fact that the derivative of a Gaussian function with respect to 25c its center position is again a Gaussian function but of a higher 26c angular momentum. Hence the same code can be used to calculate 27c derivative integrals as well as energy integrals (I am not sure 28c whether the same would be true if directly calculating the 29c derivative with respect to C as that would change the operator). 30c 31c So we expect this code to evaluate eqs. (1) and (2) and use 32c translational invariance to evaluate eq. (3) (i.e. the second 33c equality). 34c 35c Implicit real*8 (a-h,o-z) 36c Implicit integer (i-n) 37 Implicit none 38 39#include "errquit.fh" 40 41 Logical dryrun 42 Logical canAB,GenCon ! what is canAB??? 43 Integer Nint, maxW0 44 45c--> Cartesian Coordinates, Primitives & Contraction Coefficients 46 47 Double precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) 48 Double precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) 49 50c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents 51 52 Double precision Cxyz(3), zan 53 54c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals 55 56 Double precision bNAI(Nint,9) 57 58c--> Derivative Indices 59 60 Integer inder1(3,3),inder2(6,6) 61 62c--> Scratch Space. 63 64 Double precision W0(maxW0) 65 66c--> Other stuff 67 68 Integer lprod, mxd, Lp3, Lp, Nder, NPP 69 Integer NPA, NCA, NPB, NCB, NCP 70 Integer La, Lb, Li, ictrA, ictrB 71 Integer i_VA, i_VB, i_VP, i_VR, i_Ri, i_Rj, i_RS 72 Integer MaxMem, nd 73 Integer i_P, i_R0, i_top, i_pf, i_PC, i_R0C, i_left, i_right 74 Integer i_IPAIRp, i_ESp, i_ALPHAp, i_Ep 75 Integer i_IJK, i_ff 76 77 Data inder1/ 1,0,0, 78 & 0,1,0, 79 & 0,0,1 / 80 81 Data inder2/ 1,0,0,0,0,0, 82 & 0,1,0,0,0,0, 83 & 0,0,1,0,0,0, 84 & 0,0,0,1,0,0, 85 & 0,0,0,0,1,0, 86 & 0,0,0,0,0,1 / 87c 88c Compute gradient of the overlap, kinetic energy, and nuclear attraction 89c integrals for 2 shells of contracted Gaussians functions. 90c 91c The derivative indices specify the order of differentiation for each 92c coordinate: 93c 94c inder1(i,n), i=1-3: d/dRx, d/dRy, d/dRz 95c 96c inder2(i,n), i=1-3: d/dPx, d/dPy, d/dPz 97c i=4-6: d/dRx, d/dRy, d/dRz 98c 99c****************************************************************************** 100 101#if defined(INTDEBUG) 102 call hf_print_set(1) 103 call hf_print('hf1d: a shell',axyz,aprims,acoefs,npa,nca,la) 104 call hf_print('hf1d: b shell',bxyz,bprims,bcoefs,npb,ncb,lb) 105 call hf_print_set(0) 106#endif 107 MXD = 1 108 Li = 0 109 110c Determine whether general or segmented contraction is used. 111 112 NCP = NCA*NCB 113 114 GenCon = NCP.ne.1 115 116 if( GenCon )then 117 write(*,*) 'HF1D: Not yet ready for general contraction.' 118 call errquit('HF1D: Not yet ready for general contraction.', 119 & 0,UERR) 120 end if 121 122c Define the angular momentum of the overlap distribution. 123 124 Lp = La + Lb 125 126c Increment "Lp" to account for the order of differentiation. 127 128 Lp = Lp + MXD 129 130c Define the accumulated number of angular momentum functions <= Lp. 131 132 Lp3 = ((Lp+1)*(Lp+2)*(Lp+3))/6 133 134c Define the prefactor of the overlap distribution "P". 135 136c Assign pointers to scratch space. 137 138 i_ALPHAp = 1 139 i_IPAIRp = i_ALPHAp + 2*(NPA*NPB) 140 i_left = i_IPAIRp + 2*(NPA*NPB) - 1 141 142 i_ESp = (maxW0+1) - 3*(NPA*NPB) 143 i_right = i_ESp 144 145 if( i_left.ge.i_right )then 146 147 write(*,*) 'HF1D 1:Insufficient scratch space.' 148 write(*,*) ' needed ',i_left + (maxW0-(i_right-1)) 149 write(*,*) ' allocated ',maxW0 150 151 write(*,*) 'From the left ' 152 write(*,*) 'ALPHAp: ',i_ALPHAp 153 write(*,*) 'IPAIRp: ',i_IPAIRp 154 write(*,*) 'From the right ' 155 write(*,*) 'ESp : ',i_ESp 156 157 call errquit('HF1D: Insufficient scratch space.',i_right-i_left, 158 & UERR) 159 160 end if 161 162 MaxMem = 1 ! take care of compiler warnings 163 if (dryrun) then 164 MaxMem = i_left + (maxW0 - (i_right-1)) 165 NPP = NPA*NPB 166 else 167 call hfset(Axyz,Aprims,Acoefs,NPA,NCA, 168 & Bxyz,Bprims,Bcoefs,NPB,NCB, 169 & GenCon,W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ESp),NPP) 170 endif 171 172c Define the Hermite linear expansion coefficients. 173 174c Assign pointers to scratch space. 175 176 lprod = ((La+Li)+(Lb+Li)+1)*((La+Li)+1)*((Lb+Li)+1) 177 178 i_Ep = i_IPAIRp + 2*(NPA*NPB) 179 i_pf = i_Ep + 3*NPP*(MXD+1)*lprod 180 i_left = i_pf + 2*NPP - 1 181 182 if( i_left.ge.i_right )then 183 184 write(*,*) 'HF1D 2:Insufficient scratch space.' 185 write(*,*) ' needed ',i_left + (maxW0-(i_right-1)) 186 write(*,*) ' allocated ',maxW0 187 188 write(*,*) 'From the right ' 189 write(*,*) 'ALPHAp: ',i_ALPHAp 190 write(*,*) 'IPAIRp: ',i_IPAIRp 191 write(*,*) 'Ep : ',i_Ep 192 write(*,*) 'pf : ',i_pf, npp,mxd,lprod 193 write(*,*) 'From the left ' 194 write(*,*) 'ESp : ',i_ESp 195 196 call errquit('HF1D: Insufficient scratch space.',i_right-i_left, 197 & UERR) 198 199 end if 200 201 if (dryrun) then 202 MaxMem = max(MaxMem, (i_left+(maxW0-(i_right-1)))) 203 else 204 do 100 nd = 0,MXD 205 206 call hfmke(Axyz,Bxyz,W0(i_ALPHAp),W0(i_ESp),W0(i_Ep), 207 & W0(i_pf),nd,NPP,MXD,La+Li,Lb+Li) 208 209100 continue 210 endif 211 212 213 214c Compute nuclear attraction integrals, <a|V|b>. 215 216c if( NAI )then 217 218c Define the auxiliary function integrals. 219 220c Assign scratch space. 221 222 i_R0 = i_Ep + 3*NPP*(MXD+1)*lprod 223 i_IJK = i_R0 + NPP*Lp3 224 i_R0C = i_IJK + (Lp+1)**3 225 i_P = i_R0C + NPP*Lp3*1 226 i_RS = i_P + NPP*3 227 i_PC = i_RS + NPP 228 i_ff = i_PC + NPP*3 229 i_Rj = i_ff + NPP*2 230 i_top = i_Rj + NPP*(Lp+1)*Lp3 - 1 231 232 if( i_top.gt.maxW0 )then 233 234 write(*,*) 'HF1D 3:Insufficient scratch space.' 235 write(*,*) ' needed ',i_top 236 write(*,*) ' allocated ',maxW0 237 238 write(*,*) 'ALPHAp: ',i_ALPHAp 239 write(*,*) 'IPAIRp: ',i_IPAIRp 240 write(*,*) 'Ep : ',i_Ep 241 write(*,*) 'R0 : ',i_R0 242 write(*,*) 'IJK : ',i_IJK 243 write(*,*) 'R0C : ',i_R0C 244 write(*,*) 'P : ',i_P 245 write(*,*) 'RS : ',i_RS 246 write(*,*) 'PC : ',i_PC 247 write(*,*) 'ff : ',i_ff 248 write(*,*) 'Rj : ',i_Rj 249 250 call errquit('HF1D: Insufficient scratch space.',maxW0-i_top, 251 & UERR) 252 253 end if 254 255 if (dryrun) then 256 MaxMem = max(MaxMem, i_top) 257 else 258 call hf1mkr(Axyz,Bxyz,Cxyz,zan,0.0d0,1, 259 & W0(i_ALPHAp),W0(i_P),W0(i_RS),W0(i_PC),W0(i_ff), 260 & W0(i_Rj),W0(i_R0),W0(i_R0C),W0(i_IJK), 261 & NPP,Lp,Lp3,.TRUE.) 262c 263c Compute the 1-electron - nuclear attraction contributions of 264c the gradient wrt the nuclear coordinates 265c 266c call hfefi(W0(i_Ep),W0(i_R0C),W0(i_IJK),bNAI, 267c & NPP,Nint,La,Lb,Li,Lp,Lp3,1, 268c & MXD,canAB,ictra,ictrb) 269 270 endif 271c Compute compute the derivative wrt to the centers of the Gaussian functions. 272 273c Assign scratch space. 274 275 i_VP = i_IJK + (Lp+1)**3 276 i_VR = i_VP + NPP*(Nint*3) 277 i_VA = i_VR + NPP*(Nint*3) 278 i_VB = i_VA + (Nint*3) 279 i_ff = i_VB + (Nint*3) 280 i_top = i_ff + NPP*2 - 1 281 282 if( i_top.gt.maxW0 )then 283 284 write(*,*) 'HF1D 4:Insufficient scratch space.' 285 write(*,*) ' needed ',i_top 286 write(*,*) ' allocated ',maxW0 287 288 write(*,*) 'ALPHAp: ',i_ALPHAp 289 write(*,*) 'IPAIRp: ',i_IPAIRp 290 write(*,*) 'Ep : ',i_Ep 291 write(*,*) 'R0 : ',i_R0 292 write(*,*) 'IJK : ',i_IJK 293 write(*,*) 'VP : ',i_VP 294 write(*,*) 'VR : ',i_VR 295 write(*,*) 'VA : ',i_VA 296 write(*,*) 'VB : ',i_VB 297 write(*,*) 'ff : ',i_ff 298 299 call errquit('HF1D: Insufficient scratch space.',maxW0-i_top, 300 & UERR) 301 302 end if 303 304 if (dryrun) then 305 MaxMem = max(MaxMem, i_top) 306 else 307c Compute the derivatives wrt to (P,R). 308 309c Compute the derivatives of the primitive integrals. 310 311 Nder = 6 312 313 call hfdnai(W0(i_Ep),W0(i_R0),W0(i_IJK),W0(i_VP), 314 & NPP,Nint,La,Lb,Li,Lp,Lp3,1, 315 & MXD,inder2,Nder,canAB) 316 317c Transform to obtain derivatives wrt (A,B). 318 319 call hf1PRtoAB(W0(i_VP),W0(i_VR),W0(i_VA),W0(i_VB), 320 & W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ff),NPP,Nint*3, 321 & ictrA,ictrB) 322 323c Combine the two components of the derivatives of NAIs. 324c 325c (1) the derivative wrt to each of the nuclear attraction centers, 326c (2) the derivative wrt to the centers of the Gaussian functions. 327 328 if( ictrA.eq.ictrB )then 329 call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,1),1) 330 call daxpy(Nint*3,-1.D0,W0(i_VA),1,bNAI(1,7),1) 331 else 332 call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,1),1) 333 call daxpy(Nint*3,1.D0,W0(i_VB),1,bNAI(1,4),1) 334 call daxpy(Nint*3,-1.D0,W0(i_VA),1,bNAI(1,7),1) 335 call daxpy(Nint*3,-1.D0,W0(i_VB),1,bNAI(1,7),1) 336 end if 337 338 endif 339c end if 340 341c Return the maximum amount of scratch space required by a "dry run". 342 343 if( DryRun ) maxW0 = MaxMem 344c 345 end 346