1C $Id$ 2************************************************************************ 3* * 4c:tex-% this is part of the API Standard Integral routines. 5c:tex-\subsection{ecp\_integral} 6c:tex-this routine computes the scalar and spin-orbit ECP integrals. The 7c:tex-ECPs are defined by 8c:tex-\begin{eqnarray*} 9c:tex-U^{AREP} = U^{AREP}_L + \sum_{l=0}^{L-1} U^{AREP}_\ell 10c:tex-\sum_{m=-\ell}^\ell | \ell m \rangle \langle \ell m | \\ 11c:tex-U^{so} = \sum_{\ell=1}^L \frac1{\ell} \frac{2}{2\ell+1} 12C:tex-\Delta U^{so}_\ell 13c:tex-\sum_{m=-\ell}^\ell \sum_{m'=-\ell}^\ell | \ell m \rangle 14c:tex-\langle \ell m | -i\hat\ell | \ell m' \rangle \langle \ell m' | 15c:tex-\end{eqnarray*} 16c:tex-Note that the spin-orbit operator is over $-i\ell$ and that the 17c:tex-scalar product with $\sigma$ not {\bf s} should be taken. 18c:tex-Note also that the integrals come out in the order z,y,x. 19c:tex- 20c:tex- 21c:tex- 22c:tex-{\it Syntax:} 23c:tex-\begin{verbatim} 24 subroutine ecp_integral ( 25 & xyz_A,zeta_A,coef_A,n_prim_A,n_cont_A,l_A,i_c_A, 26 & xyz_B,zeta_B,coef_B,n_prim_B,n_cont_B,l_B,i_c_B, 27 & xyz_C,zeta_C,coef_C,n_prim_C,n_colc_C, 28 & ind_z,ind_c,n_zeta_C,n_coef_C, 29 & l_C,i_cent_C,n_C,l_ecp_max, 30 & sphcart,csco,lcsco, 31 & ecp_ints,n_int,n_blk, 32 & DryRun,scr,lscr,ibug) 33c:tex-\end{verbatim} 34* * 35* Calculate scalar and spin-orbit ecp integrals. * 36* * 37* Argument (status) - description * 38* * 39* xyz_A (inp) - coordinates of centre A * 40* zeta_A (inp) - exponents of primitive gaussians on centre A * 41* coef_A (inp) - contraction coefficients on centre A * 42* n_prim_A (inp) - number of primitive gaussians on centre A * 43* n_cont_A (inp) - number of contracted functions on centre A * 44* l_A (inp) - angular momentum of functions on centre A * 45* i_c_A (inp) index of centre A in centre list * 46* xyz_B (inp) - coordinates of centre B * 47* zeta_B (inp) - exponents of primitive gaussians on centre B * 48* coef_B (inp) - contraction coefficients on centre B * 49* n_prim_B (inp) - number of primitive gaussians on centre B * 50* n_cont_B (inp) - number of contracted functions on centre B * 51* l_B (inp) - angular momentum of functions on centre B * 52* i_c_B (inp) index of centre B in centre list * 53* xyz_C (inp) - coordinates of ECP centres C * 54* zeta_C (inp) - array of exponents of primitive gaussians on all * 55* centres C. These are stored in an array of single * 56* dimension, i.e. packed. * 57* coef_C (inp) - array of contraction coefficients on all centres C * 58* n_prim_C (inp) - array of number of primitive gaussians for each * 59* power of r, l value and ECP centre. The highest * 60* l value is for the local part, thus the second * 61* dimension is l_ecp_max+2 (or 0:l_ecp_max+1) * 62* n_colc_C (inp) - array of number of coefficients for each l value * 63* and ECP centre. This is n_prim_C summed over the * 64* first dimension. * 65* ind_z (inp) - array of addresses of first exponent for each l * 66* value and ECP centre. * 67* ind_c (inp) - array of addresses of first coefficient for each l * 68* value and ECP centre. * 69* n_zeta_C (inp) - total number of ECP exponents. * 70* n_coef_C (inp) - total number of ECP coefficients. * 71* l_C (inp) - maximum angular momentum of projectors on centres C * 72* i_cent_C (inp) - indices of ECP centres C * 73* n_C (inp) - number of ECP centres C * 74* l_ecp_max (inp) - maximum angular momentum of any projector on any * 75* ECP centre. * 76* sphcart - 0 for cartesian integrals, 1 for spherical integrals * 77* ecp_ints (out) - integrals over ECPs * 78* n_int (inp) - number of ECP integrals. Should be equal to * 79* NCA*NCB*[(La+1)*(La+2)/2]*[(Lb+1)*(Lb+2)/2] * 80* n_blk (inp) - 1 for scalar only, 3 for s-o only, 4 for both * 81* s-o integrals come out in the order z,y,x. * 82* DryRun (inp) - logical for dry run. If true, routine only returns * 83* maximum scratch space needed, if false, integrals * 84* are returned. * 85* scr (scr) - scratch array for work space * 86* lscr (i/o) - length of scratch array. Value returned if DryRun is * 87* true, used as dimension if false. * 88* ibug - debug flag. 0 for no debug, 1 for address printing, 2 for * 89* array printing, 3 for both. * 90* * 91* Written by K. G. Dyall * 92* * 93************************************************************************ 94 implicit none 95#include "stdio.fh" 96#include "ecp_consts.fh" 97#include "util.fh" 98#include "errquit.fh" 99 integer i,j,l,n, 100 & n_prim_A,n_cont_A,l_A,i_c_A, 101 & n_prim_B,n_cont_B,l_B,i_c_B, 102 & n_zeta_C,n_coef_C,n_C,i_c_C,l_ecp_max, 103 & n_int,n_blk,sphcart,lscr,ibug 104 integer n_prim_C(0:4,-1:l_ecp_max,n_C*2), 105 & n_colc_C(-1:l_ecp_max,n_C*2), 106 & ind_z(-1:l_ecp_max,n_C*2),ind_c(-1:l_ecp_max,n_C*2), 107 & l_C(n_C),i_cent_C(n_C) 108 integer i_scr,i_coef,i_zeta,memscr,mem, 109 & max_type0,max_type1,max_type2,n_intt,n_all_a,n_all_b, 110 & n_cart_a,n_cart_b,n_cart_ab,n_cont_ab 111 integer lcsco 112 logical DryRun,debug_gen,debug_addresses,debug_arrays 113 double precision 114 & xyz_A(3),zeta_A(n_prim_A),coef_A(n_prim_A,n_cont_A), 115 & xyz_B(3),zeta_B(n_prim_B),coef_B(n_prim_B,n_cont_B), 116 & xyz_C(3,n_C),zeta_C(n_zeta_C),coef_C(n_zeta_C), 117 & scr(lscr),ecp_ints(n_int*n_blk) 118 double precision 119 & X_AC,Y_AC,Z_AC,R_AC,X_BC,Y_BC,Z_BC,R_BC, 120 & tol 121 double precision csco(lcsco) 122 data tol/1.0d-16/ 123* 124 if ((n_blk .ne. 1) .and. (n_blk .ne. 3) .and. (n_blk .ne. 4)) 125 & call errquit('Illegal value of n_blk in ecp_integral',99, 126 & INT_ERR) 127* 128 debug_gen = ibug .gt. 0 129 debug_addresses = mod(ibug,2) .eq. 1 130 debug_arrays = mod(ibug,10)/2 .eq. 1 131 if (debug_gen) write (LuOut,*) ibug 132* 133 if (debug_gen) write (LuOut,'(//A,/)') 'Entering ecp_integral ...' 134 if (debug_addresses) then 135 write(LuOut,*)' lscr in ecp_integral:',lscr 136 write (LuOut,*) 'n_prim_A,n_cont_A,l_A',n_prim_A,n_cont_A,l_A 137 write (LuOut,*) 'n_prim_B,n_cont_B,l_B',n_prim_B,n_cont_B,l_B 138 write (LuOut,*) 'l_ecp_max,n_c',l_ecp_max,n_c 139 end if 140 n_cart_a = (l_a+1)*(l_a+2)/2 141 n_cart_b = (l_b+1)*(l_b+2)/2 142 n_cart_ab = n_cart_a*n_cart_b 143 n_cont_ab = n_cont_a*n_cont_b 144 n_all_a = n_cart_a*n_cont_a 145 n_all_b = n_cart_b*n_cont_b 146 n_intt = n_cart_ab*n_cont_ab 147 if (debug_addresses) write (LuOut,*) 148 & 'n_cart_a,n_cart_b,n_cart_ab,n_cont_ab', 149 & n_cart_a,n_cart_b,n_cart_ab,n_cont_ab 150 if (sphcart .ne. 0) call errquit( 151 & 'Do your own spherical transformation, lazy bum!',99, 152 & BASIS_ERR) 153 if (n_int .lt. n_intt) then 154 write (LuOut,*) 'n_int =',n_int 155 write (LuOut,*) 'n_intt =',n_intt 156 call errquit ( 157 & 'Mismatch of integral count in ecp_integrals',99, 158 & BASIS_ERR) 159 end if 160 i_scr = n_intt*n_blk+1 161 memscr = lscr-i_scr+1 162 if (DryRun) then 163 max_type0 = 0 164 max_type1 = 0 165 max_type2 = 0 166 else 167 call dcopy (n_intt*n_blk,zero,0,ecp_ints,1) 168 if (debug_addresses) write (LuOut,*) 'memscr =',memscr 169 if (memscr .lt. 0) call errquit ( 170 & 'Insufficient scratch memory in ecp_integral',99, MEM_ERR) 171 end if 172 if (debug_addresses) write (LuOut,*) 'i_scr',i_scr 173 if (debug_arrays) then 174 call ecp_matpr(xyz_A,1,3,1,1,1,3,1,1,'Centre A coordinates', 175 & 'E',78,4) 176 call ecp_matpr(xyz_B,1,3,1,1,1,3,1,1,'Centre B coordinates', 177 & 'E',78,4) 178 call ecp_matpr(xyz_C,1,3,1,n_C,1,3,1,n_C,'ECP coordinate array', 179 & 'E',78,4) 180 call ecp_matpr(zeta_A,1,n_prim_A,1,1,1,n_prim_A,1,1, 181 & 'A exponents','E',78,4) 182 call ecp_matpr(coef_A,1,n_prim_A,1,n_cont_A,1,n_prim_A, 183 & 1,n_cont_A,'A coefficients','E',78,4) 184 call ecp_matpr(zeta_B,1,n_prim_B,1,1,1,n_prim_B,1,1, 185 & 'B exponents','E',78,4) 186 call ecp_matpr(coef_B,1,n_prim_B,1,n_cont_B,1,n_prim_B, 187 & 1,n_cont_B,'B coefficients','E',78,4) 188 call ecp_matpr(zeta_C,1,n_zeta_C,1,1,1,n_zeta_C,1,1, 189 & 'ECP exponents','E',78,4) 190 call ecp_matpr(coef_C,1,n_zeta_C,1,1,1,n_zeta_C,1,1, 191 & 'ECP coefficients','E',78,4) 192 end if 193 if (debug_gen) write (LuOut,*) 'Number of ECP centers =',n_C 194* 195* Loop over ECP centres 196* 197 do i = 1,n_C 198 l = l_C(i) 199 i_c_C = i_cent_C(i) 200 if (debug_gen) write (LuOut,*) 'ECP center',i 201 if (debug_gen) write (LuOut,*) ' Maximum angular momentum',l 202* 203* Set up relative cartesian coordinates 204* 205 X_AC = xyz_C(1,i)-xyz_A(1) 206 Y_AC = xyz_C(2,i)-xyz_A(2) 207 Z_AC = xyz_C(3,i)-xyz_A(3) 208 R_AC = sqrt(X_AC**2+Y_AC**2+Z_AC**2) 209 X_BC = xyz_C(1,i)-xyz_B(1) 210 Y_BC = xyz_C(2,i)-xyz_B(2) 211 Z_BC = xyz_C(3,i)-xyz_B(3) 212 R_BC = sqrt(X_BC**2+Y_BC**2+Z_BC**2) 213 if (debug_arrays) then 214 write (LuOut,'(3x,A,3F10.6)') 'Relative coords of center A:', 215 & X_AC,Y_AC,Z_AC 216 write (LuOut,'(3x,A,3F10.6)') 'Relative coords of center B:', 217 & X_BC,Y_BC,Z_BC 218 write (LuOut,'(3x,A,3F10.6)') 'Distance to centers A and B:', 219 & R_AC,R_BC 220 end if 221* 222* Evaluate integrals 223* 224 if (debug_gen) write (LuOut,*) 'starting integrals' 225C if (debug_arrays) then 226C write (LuOut,*) 'Parameters of ECP' 227C call ecp_matpi (n_prim_C,0,4,-1,l_ecp_max,0,2,0,l_C, 228C & 'n_prim_C',81,15) 229C call ecp_matpi (ind_C,-1,l_ecp_max,1,2*n_C,0,l_C,1,2*n_C, 230C & 'ind_C',81,15) 231C call ecp_matpi (n_colc_C,-1,l_ecp_max,1,2*n_C,0,l_C,1,2*n_C, 232C & 'n_colc_C',81,15) 233C call ecp_matpr (zeta_C,1,n_zeta_C,1,1,i_zeta,i_zeta+n-1,1,1, 234C & 'zeta_C','E',81,5) 235C call ecp_matpr (coef_C,1,n_zeta_C,1,1,i_zeta,i_zeta+n-1,1,1, 236C & 'coef_C','E',81,5) 237C end if 238 if (.not.DryRun) call dcopy (n_intt*n_blk,zero,0,scr,1) 239 i_zeta = ind_z(-1,i) 240 i_coef = ind_c(-1,i) 241 if ((i_c_A .eq. i_c_C) .and. (i_c_B .eq. i_c_C)) then 242* 243* One-centre case, A = B = C 244* 245 if (debug_gen) write (LuOut,*) 'One-centre case' 246 if ((n_colc_C(-1,i) .gt. 0)) then 247 if ((n_blk .ne. 3)) 248 & call ecp_local0 (mem,DryRun, 249 & l_B,n_prim_B,n_cont_B,coef_B,zeta_B,n_cart_b, 250 & l_A,n_prim_A,n_cont_A,coef_A,zeta_A,n_cart_a, 251 & n_prim_C(0,-1,i),n_colc_C(-1,i), 252 & zeta_c(i_zeta),coef_c(i_coef),0,4, 253 & tol,sphcart,scr(i_scr),memscr, 254 & csco,lcsco, 255 & scr,ibug/10) 256 if (debug_arrays .and. .not.DryRun) 257 & call ecp_matpr (scr,1,n_all_b,1,n_all_a,1,n_all_b, 258 & 1,n_all_a,'Local ECP integrals','E',78,4) 259 if (DryRun) max_type0 = max(mem,max_type0) 260 end if 261 call ecp_int0 (mem,DryRun, 262 & l_B,n_prim_B,n_cont_B,coef_B,zeta_B,n_cart_b, 263 & l_A,n_prim_A,n_cont_A,coef_A,zeta_A,n_cart_a, 264 & l_C(i),n_prim_C(0,-1,i),n_colc_C(-1,i),ind_z(-1,i), 265 & ind_c(-1,i),n_zeta_C,n_coef_C,l_ecp_max,n_C, 266 & zeta_C,coef_C,0,4,tol,sphcart, 267 & scr(i_scr),memscr,csco,lcsco, 268 & scr,n_blk,ibug/10) 269 if (DryRun) max_type0 = max(mem,max_type0) 270 else if (i_c_A .eq. i_c_C) then 271* 272* Two-centre case with A = C 273* 274 if (debug_gen) write (LuOut,*) 'Two-centre case A=C' 275 if ((n_colc_C(-1,i) .gt. 0)) then 276 if ((n_blk .ne. 3)) 277 & call ecp_local1 (mem,DryRun, 278 & R_BC,X_BC,Y_BC,Z_BC,l_B,n_prim_B,n_cont_B,coef_B,zeta_B, 279 & R_AC,X_AC,Y_AC,Z_AC,l_A,n_prim_A,n_cont_A,coef_A,zeta_A, 280 & n_prim_C(0,-1,i),n_colc_C(-1,i), 281 & zeta_c(i_zeta),coef_c(i_coef),0,4, 282 & tol,sphcart,scr(i_scr),memscr, 283 & csco,lcsco, 284 & scr,n_intt,.false.,ibug/10) 285 if (debug_arrays .and. .not.DryRun) 286 & call ecp_matpr (scr,1,n_all_b,1,n_all_a,1,n_all_b, 287 & 1,n_all_a,'Local ECP integrals','E',78,4) 288 if (DryRun) max_type1 = max(mem,max_type1) 289 end if 290 call ecp_int1 (mem,DryRun, 291 & R_BC,X_BC,Y_BC,Z_BC,l_B,n_prim_B,n_cont_B,coef_B,zeta_B, 292 & R_AC,X_AC,Y_AC,Z_AC,l_A,n_prim_A,n_cont_A,coef_A,zeta_A, 293 & l_C(i),n_prim_C(0,-1,i),n_colc_C(-1,i),ind_z(-1,i), 294 & ind_c(-1,i),n_zeta_C,n_coef_C,l_ecp_max,n_C, 295 & zeta_C,coef_C,0,4,tol,sphcart,scr(i_scr),memscr, 296 & csco,lcsco, 297 & scr,n_intt,n_blk,.false.,ibug/10) 298 if (DryRun) max_type1 = max(mem,max_type1) 299 else if (i_c_B .eq. i_c_C) then 300* 301* Two-centre case with B = C 302* 303 if (debug_gen) write (LuOut,*) 'Two-centre case B=C' 304 if ((n_colc_C(-1,i) .gt. 0)) then 305 if ((n_blk .ne. 3)) 306 & call ecp_local1 (mem,DryRun, 307 & R_AC,X_AC,Y_AC,Z_AC,l_A,n_prim_A,n_cont_A,coef_A,zeta_A, 308 & R_BC,X_BC,Y_BC,Z_BC,l_B,n_prim_B,n_cont_B,coef_B,zeta_B, 309 & n_prim_C(0,-1,i),n_colc_C(-1,i), 310 & zeta_c(i_zeta),coef_c(i_coef),0,4, 311 & tol,sphcart,scr(i_scr),memscr, 312 & csco,lcsco, 313 & scr,n_intt,.true.,ibug/10) 314 if (debug_arrays .and..not.DryRun) 315 & call ecp_matpr (scr,1,n_all_b,1,n_all_a,1,n_all_b, 316 & 1,n_all_a,'Local ECP integrals','E',78,4) 317 if (DryRun) max_type1 = max(mem,max_type1) 318 end if 319 call ecp_int1 (mem,DryRun, 320 & R_AC,X_AC,Y_AC,Z_AC,l_A,n_prim_A,n_cont_A,coef_A,zeta_A, 321 & R_BC,X_BC,Y_BC,Z_BC,l_B,n_prim_B,n_cont_B,coef_B,zeta_B, 322 & l_C(i),n_prim_C(0,-1,i),n_colc_C(-1,i),ind_z(-1,i), 323 & ind_c(-1,i),n_zeta_C,n_coef_C,l_ecp_max,n_C, 324 & zeta_C,coef_C,0,4,tol,sphcart,scr(i_scr),memscr, 325 & csco,lcsco, 326 & scr,n_intt,n_blk,.true.,ibug/10) 327 if (DryRun) max_type1 = max(mem,max_type1) 328 else 329* 330* Three-centre case 331* 332 if (debug_gen) write (LuOut,*) 'Three-centre case' 333 if ((n_colc_C(-1,i) .gt. 0))then 334 if ((n_blk .ne. 3)) 335 & call ecp_local2 (mem,DryRun, 336 & R_BC,X_BC,Y_BC,Z_BC,l_B,n_prim_B,n_cont_B,coef_B,zeta_B, 337 & R_AC,X_AC,Y_AC,Z_AC,l_A,n_prim_A,n_cont_A,coef_A,zeta_A, 338 & n_prim_C(0,-1,i),n_colc_C(-1,i), 339 & zeta_C(i_zeta),coef_C(i_coef),0,4, 340 & tol,sphcart,scr(i_scr),memscr, 341 & csco,lcsco, 342 & scr,n_intt,ibug/10) 343 if (debug_arrays .and..not.DryRun) 344 & call ecp_matpr (scr,1,n_all_b,1,n_all_a,1,n_all_b, 345 & 1,n_all_a,'Local ECP integrals','E',78,4) 346 if (DryRun) max_type2 = max(mem,max_type2) 347 end if 348 call ecp_int2 (mem,DryRun, 349 & R_BC,X_BC,Y_BC,Z_BC,l_B,n_prim_B,n_cont_B,coef_B,zeta_B, 350 & R_AC,X_AC,Y_AC,Z_AC,l_A,n_prim_A,n_cont_A,coef_A,zeta_A, 351 & l_c(i),n_prim_C(0,-1,i),n_colc_C(-1,i),ind_z(-1,i), 352 & ind_c(-1,i),n_zeta_C,n_coef_C,l_ecp_max,n_C, 353 & zeta_C,coef_C,0,4,tol,sphcart,scr(i_scr),memscr, 354 & csco,lcsco, 355 & scr,n_intt,n_blk,ibug/10) 356 if (DryRun) max_type2 = max(mem,max_type2) 357 end if 358* 359* Add integrals from this centre to the total. 360* 361 if (.not.DryRun) then 362 call daxpy (n_intt*n_blk,one,scr,1,ecp_ints,1) 363 if (debug_arrays) then 364 n = 1 365 do j = 1,n_blk 366 call ecp_matpr (scr(n),1,n_all_b,1,n_all_a,1,n_all_b, 367 & 1,n_all_a,'Non-local ECP integrals','E',78,4) 368 n = n+n_intt 369 end do 370 end if 371 end if 372 end do 373 if (DryRun) then 374 lscr = max(max_type0,max_type1,max_type2)+n_intt*n_blk 375 if (debug_addresses) then 376 write (LuOut,*) 'max_type0',max_type0 377 write (LuOut,*) 'max_type1',max_type1 378 write (LuOut,*) 'max_type2',max_type2 379 write (LuOut,*) 'n_intt',n_intt 380 write (LuOut,*) 'lscr',lscr 381 end if 382 else if (debug_arrays) then 383 n = 1 384 do j = 1,n_blk 385 call ecp_matpr (ecp_ints(n),1,n_all_b,1,n_all_a, 386 & 1,n_all_b,1,n_all_a,'Final ECP integrals','E',78,4) 387 n = n+n_intt 388 end do 389 end if 390 if (debug_gen) write (LuOut,*) 'Exiting ecp_integral' 391* 392 return 393 end 394