1C $Id$ 2************************************************************************ 3c:tex-\subsection{rel\_oneld} 4c:tex-This routine generates the modified one-electron gradient integrals 5c:tex-for a relativistic basis set. These are the gradients of the modified 6c:tex-kinetic energy, the modified potential energy and the modified overlap, 7c:tex-\begin{eqnarray} 8c:tex-&& \tilde{T}_{ab} = T_{ab}^{LS} + T_{ab}^{SL} - T_{ab}^{SS} 9c:tex- \nonumber \\ 10c:tex-&& \tilde{V}^{sf}_{ab} = V_{ab}^{LL} + {{\alpha^2}\over{4}} 11c:tex- \nabla_A\cdot\nabla_B V_{ab}^{SS} 12c:tex-&& \tilde{V}^{so}_{ab} = V_{ab}^{LL} + {{\alpha^2}\over{4}} 13c:tex- \nabla_A\times\nabla_B V_{ab}^{SS} 14c:tex- \nonumber \\ 15c:tex-&& \tilde{S}_{ab} = S_{ab}^{LL} + {{\alpha^2}\over{2}} T_{ab}^{SS} 16c:tex- \nonumber 17c:tex-\end{eqnarray} 18c:tex- 19c:tex-\noindent Author: K. G. Dyall 20c:tex- 21c:tex-{\it Syntax:} 22c:tex-\begin{verbatim} 23 subroutine rel_oneld_cosmo ( 24 & Axyz,zeta_A,coefL_A,coefS_A,n_prim_A,n_cont_A,l_A,ictr_A, 25 & Bxyz,zeta_B,coefL_B,coefS_B,n_prim_B,n_cont_B,l_B,ictr_B, 26 & Cxyz,zan,V,lstv,canAB, 27 & do_nw,do_hnd,nonrel,DryRun,scr,lscr,ibug,ntyp) 28c:tex-\end{verbatim} 29 implicit none 30#include "stdio.fh" 31#include "rel_consts.fh" 32#include "errquit.fh" 33* 34c:tex-{\it Argument list:} 35c:tex-\begin{verbatim} 36 integer n_prim_A ! [input] num. prims in shell A 37 integer n_cont_A ! [input] num general conts in shell A 38 integer l_A ! [input] angular momentum of shell A 39 integer ictr_A ! [input] lexical atom index for shell A 40 integer n_prim_B ! [input] num. prims in shell B 41 integer n_cont_B ! [input] num general conts in shell B 42 integer l_B ! [input] angular momentum of shell B 43 integer ictr_B ! [input] lexical atom index for shell B 44 integer lscr ! [input] size of scratch array 45 integer lstv ! [input] size of any integral buffer 46 integer ibug ! [input] debug variable 47 integer ntyp ! [input] potential energy integral type 48 double precision Axyz(3) ! [input] position of center A 49 double precision zeta_A(n_prim_A) ! [input] exponents of shell A 50 double precision coefL_A(n_prim_A,n_cont_A) ! [input] A large coeffs 51 double precision coefS_A(n_prim_A,n_cont_A) ! [input] A small coeffs 52 double precision Bxyz(3) ! [input] position of center B 53 double precision zeta_B(n_prim_B) ! [input] exponents of shell B 54 double precision coefL_B(n_prim_B,n_cont_B) ! [input] B large coeffs 55 double precision coefS_B(n_prim_B,n_cont_B) ! [input] B small coeffs 56 double precision Cxyz(3,1) ! [input] all atom positions 57 double precision zan(1) ! [input] charges on all atoms 58 double precision scr(lscr) ! [scratch] scratch buffers 59 double precision V(lstv*3*3,ntyp) ! [output] potential integrals 60 logical canAB ! [input] compute only canonical ints (false only) 61 logical do_nw ! [input] can do NW integrals 62 logical do_hnd ! [input] can do HONDO integrals 63 logical nonrel ! [input] true if either centre is nonrelativistic 64 logical DryRun ! [input] true means only compute required memory 65c:tex-\end{verbatim} 66c:tex-See rel_pot for a description of the allowed values of ibug and ntyp 67c:tex-Note that in the current version of this routine, the call to rel_pot 68c:tex-uses a dummy ntyp=1. It is kept in the input so that in future, the 69c:tex-spin-orbit integrals can also be obtained with a call to this routine. 70c:tex- 71c:tex-{\it Subroutines called:} hf1d, rel\_pot, daxpy, dcopy 72* 73 integer n_cart_a ! cartesian components of shell A 74 integer n_cart_b ! cartesian components of shell B 75 integer n_cart_ab ! n_cart_a*n_cart_b 76 integer n_cont_ab ! n_cont_a*n_cont_b 77 integer n_all_b ! n_cart_b*n_cont_b 78 integer n_all_a ! n_cart_a*n_cont_a 79 integer n_ab ! number of integrals 80 integer n_ab6 ! n_ab*6, number of gradient integrals for T and S 81 integer n_ab3at ! n_ab*3*nat, number of gradient integrals for V 82 integer n_cartp_a ! cartesian components for l_A+1 83 integer n_cartp_b ! cartesian components for l_B+1 84 integer n_cartm_a ! cartesian components for l_A-1 85 integer n_cartm_b ! cartesian components for l_B-1 86 integer n_intpp ! number of integrals for l_A+1,l_B+1 87 integer n_intpm ! number of integrals for l_A-1,l_B+1 88 integer n_intmp ! number of integrals for l_A+1,l_B-1 89 integer n_intmm ! number of integrals for l_A-1,l_B-1 90 integer i_xca ! address in scr of exp*coef for shell A 91 integer i_xcb ! address in scr of exp*coef for shell B 92 integer i_pp ! address in scr of integrals for l_A+1,l_B+1 93 integer i_pm ! address in scr of integrals for l_A-1,l_B+1 94 integer i_mp ! address in scr of integrals for l_A+1,l_B-1 95 integer i_mm ! address in scr of integrals for l_A-1,l_B-1 96 integer i_scr ! address of free space in scr 97 integer memscr ! free space in scr 98 integer max_mem ! maximum memory used 99 integer i,j,k,l,m ! loop indices etc. 100 double precision one ! Obvious! 101 parameter (one = 1.0D0) 102* 103 integer n_allp_b ! n_cartp_b*n_cont_b 104 integer n_allp_a ! n_cartp_a*n_cont_a 105 integer n_allm_b ! n_cartm_b*n_cont_b 106 integer n_allm_a ! n_cartm_a*n_cont_a 107* 108 logical debug_gen ! do general debug printing 109 logical debug_addresses ! do address debug printing 110 logical debug_arrays ! do array debug printing 111 logical doS ! compute overlap (True/False) 112 logical doT ! compute kinetic (True/False) 113 logical doV ! compute potential (True/False) 114 logical doVtil ! compute potential (True/False) 115 character*12 pot_type(4) ! potential type labels 116 character*1 xyz(3) ! coordinate labels 117 data pot_type 118 & /' Scalar','z spin-orbit','y spin-orbit','x spin-orbit'/ 119 data xyz/'x','y','z'/ 120 data doVtil / .true. / 121* 122 debug_gen = ibug .gt. 0 123 debug_addresses = mod(ibug,2) .eq. 1 124 debug_arrays = (mod(ibug,10)/2 .eq. 1) .and. .not.DryRun 125 max_mem = 0 126* 127 if (.not.(do_nw .or. do_hnd)) call errquit 128 & ('rel_oneld: can''t do NW or HONDO integrals',99, INT_ERR) 129* 130 if (debug_gen) then 131 write (LuOut,*) 'Beginning rel_oneld_cosmo' 132 write (LuOut,*) 'l_A',l_A 133 write (LuOut,*) 'n_prim_A',n_prim_A 134 write (LuOut,*) 'n_cont_A',n_cont_A 135 write (LuOut,*) 'ictr_A',ictr_A 136 write (LuOut,*) 'l_B',l_B 137 write (LuOut,*) 'n_prim_B',n_prim_B 138 write (LuOut,*) 'n_cont_B',n_cont_B 139 write (LuOut,*) 'ictr_B',ictr_B 140c if (doStil) write (LuOut,*) 'Doing overlaps' 141c if (doTtil) write (LuOut,*) 'Doing kinetic energy' 142 if (doVtil) write (LuOut,*) 'Doing potential energy' 143 call util_flush(LuOut) 144 end if 145* 146 n_cart_a = (l_a+1)*(l_a+2)/2 147 n_cart_b = (l_b+1)*(l_b+2)/2 148 n_cart_ab = n_cart_a*n_cart_b 149 n_cont_ab = n_cont_a*n_cont_b 150 n_all_a = n_cart_a*n_cont_a 151 n_all_b = n_cart_b*n_cont_b 152 n_ab = n_cart_ab*n_cont_ab 153 n_ab6 = n_ab*6 154 n_ab3at = n_ab*3*3 155 if (lstv .lt. n_ab .and. .not.DryRun) call errquit ( 156 & 'Integral buffer length too small in rel_oneld',99, 157 & MEM_ERR) 158 if (debug_addresses) then 159 write (LuOut,*) 'n_cart_a',n_cart_a 160 write (LuOut,*) 'n_cart_b',n_cart_b 161 write (LuOut,*) 'n_cart_ab',n_cart_ab 162 write (LuOut,*) 'n_cont_ab',n_cont_ab 163 write (LuOut,*) 'n_all_a',n_all_a 164 write (LuOut,*) 'n_all_b',n_all_b 165 write (LuOut,*) 'n_ab',n_ab 166 call util_flush(LuOut) 167 end if 168 if (debug_arrays) then 169 call ecp_matpr (coefL_A,1,n_prim_a,1,n_cont_a, 170 & 1,n_prim_a,1,n_cont_a,'L coef A','E',120,6) 171 call ecp_matpr (coefS_A,1,n_prim_a,1,n_cont_a, 172 & 1,n_prim_a,1,n_cont_a,'S coef A','E',120,6) 173 call ecp_matpr (coefL_B,1,n_prim_b,1,n_cont_b, 174 & 1,n_prim_b,1,n_cont_b,'L coef B','E',120,6) 175 call ecp_matpr (coefS_B,1,n_prim_b,1,n_cont_b, 176 & 1,n_prim_b,1,n_cont_b,'S coef B','E',120,6) 177 call util_flush(LuOut) 178 end if 179* 180* Calculate large component overlap and nuclear attraction integrals 181* 182c doS = doStil 183c doT = doTtil.and.nonrel 184c doV = doVtil 185 memscr = lscr 186 if (do_nw) then 187 call hf1d_cosmo( 188 & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A,ictr_A, 189 & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B,ictr_B, 190 & Cxyz,zan,V,n_ab,canAB, 191 & DryRun,scr,memscr) 192 else 193 call errquit("rel_oneld_cosmo: need NW",0,CAPMIS_ERR) 194c call hnd_stvintd( 195c & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A,ictr_A, 196c & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B,ictr_B, 197c & Cxyz,zan,nat,S,T,V,n_ab,doS,doT,doV,scr,memscr) 198 end if 199 if (debug_arrays) then 200 i = 1 201 do j = 1,3 202 write (LuOut,'(2A)') xyz(j),' component of derivatives' 203c if (doS) call ecp_matpr(S(i),1,n_all_b,1,n_all_a, 204c & 1,n_all_b,1,n_all_a,'LL overlap','E',120,6) 205c if (doT) call ecp_matpr(T(i),1,n_all_b,1,n_all_a, 206c & 1,n_all_b,1,n_all_a,'LL kinetic','E',120,6) 207 if (doV) call ecp_matpr(V(i,1),1,n_all_b,1,n_all_a, 208 & 1,n_all_b,1,n_all_a,'LL potential','E',120,6) 209 i = i+n_ab 210 call util_flush(LuOut) 211 end do 212 end if 213 if (DryRun) max_mem = max(max_mem,memscr) 214 if (debug_gen) write (LuOut,*) 'Large component done' 215 if (nonrel) return 216* 217* Calculate kinetic energy integrals, correction to overlaps 218* 219c if (doTtil) then 220c doS = .false. 221c doT = .true. 222c doV = .false. 223c memscr = lscr-n_ab6 224c if (do_nw) then 225c call hf1d( 226c & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A,ictr_A, 227c & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B,ictr_B, 228c & Cxyz,zan,exinv,nat,S,T,V,n_ab,doS,doT,doV,canAB, 229c & DryRun,scr(n_ab6+1),memscr) 230c else 231c call hnd_stvintd( 232c & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A,ictr_A, 233c & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B,ictr_B, 234c & Cxyz,zan,nat,S,T,V,n_ab,doS,doT,doV, 235c & scr(n_ab6+1),memscr) 236c end if 237c if (DryRun) max_mem = max(max_mem,memscr) 238c if (debug_arrays) then 239c i = 1 240c do j = 1,3 241c write (LuOut,'(1x,2A)') xyz(j),' component of derivatives' 242c call ecp_matpr(T(i),1,n_all_b,1,n_all_a, 243c & 1,n_all_b,1,n_all_a,'LS kinetic','E',120,6) 244c i = i+n_ab 245c call util_flush(LuOut) 246c end do 247c end if 248c memscr = lscr-n_ab6 249c if (do_nw) then 250c call hf1d( 251c & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A,ictr_A, 252c & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B,ictr_B, 253c & Cxyz,zan,exinv,nat,S,scr,V,n_ab,doS,doT,doV,canAB, 254c & DryRun,scr(n_ab6+1),memscr) 255c else 256c call hnd_stvintd( 257c & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A,ictr_A, 258c & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B,ictr_B, 259c & Cxyz,zan,nat,S,scr,V,n_ab,doS,doT,doV, 260c & scr(n_ab6+1),memscr) 261c end if 262c if (debug_arrays) then 263c i = 1 264c do j = 1,3 265c write (LuOut,'(1x,2A)') xyz(j),' component of derivatives' 266c call ecp_matpr(scr(i),1,n_all_b,1,n_all_a, 267c & 1,n_all_b,1,n_all_a,'SL kinetic','E',120,6) 268c i = i+n_ab 269c call util_flush(LuOut) 270c end do 271c end if 272c if (DryRun) then 273c max_mem = max(max_mem,memscr+n_ab6) 274c else 275c call daxpy (n_ab6,one,scr,1,T,1) 276c end if 277c if (debug_arrays) then 278c i = 1 279c do j = 1,3 280c write (LuOut,'(1x,2A)') xyz(j),' component of derivatives' 281c call ecp_matpr(T(i),1,n_all_b,1,n_all_a, 282c & 1,n_all_b,1,n_all_a,'LS+SL kinetic','E',120,6) 283c i = i+n_ab 284c end do 285c call util_flush(LuOut) 286c end if 287c end if 288c if (doStil .or. doTtil) then 289c doS = .false. 290c doT = .true. 291c doV = .false. 292c memscr = lscr-n_ab6 293c if (do_nw) then 294c call hf1d( 295c & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A,ictr_A, 296c & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B,ictr_B, 297c & Cxyz,zan,exinv,nat,S,scr,V,n_ab,doS,doT,doV,canAB, 298c & DryRun,scr(n_ab6+1),memscr) 299c else 300c call hnd_stvintd( 301c & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A,ictr_A, 302c & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B,ictr_B, 303c & Cxyz,zan,nat,S,scr,V,n_ab,doS,doT,doV, 304c & scr(n_ab6+1),memscr) 305c end if 306c if (DryRun) then 307c max_mem = max(max_mem,memscr+n_ab6) 308c else 309c if (debug_arrays) then 310c i = 1 311c do j = 1,3 312c write (LuOut,'(1x,2A)') xyz(j),' component of derivatives' 313c call ecp_matpr(scr(i),1,n_all_b,1,n_all_a, 314c & 1,n_all_b,1,n_all_a,'SS kinetic','E',120,6) 315c i = i+n_ab 316c end do 317c call util_flush(LuOut) 318c end if 319c if (doTtil) call daxpy (n_ab6,-one,scr,1,T,1) 320c if (doStil) call daxpy (n_ab6,halsq,scr,1,S,1) 321c if (debug_arrays) then 322c i = 1 323c do j = 1,3 324c write (LuOut,'(1x,2A)') xyz(j),' component of derivatives' 325c if (doTtil) call ecp_matpr(T(i),1,n_all_b,1,n_all_a, 326c & 1,n_all_b,1,n_all_a,'LS+SL-SS kinetic','E',120,6) 327c if (doStil) call ecp_matpr(S(i),1,n_all_b,1,n_all_a, 328c & 1,n_all_b,1,n_all_a,'LL+SS overlap','E',120,6) 329c i = i+n_ab 330c call util_flush(LuOut) 331c end do 332c end if 333c end if 334c if (debug_gen) write (LuOut,*) 'KE & overlap done' 335c end if 336 if (.not.doVtil) return 337* 338* Generate small component potential arrays 339* 340* 341* Set up pointers to scratch space for coefficients multiplied by 342* exponents and for integrals with shifted l values 343* 344 n_cartp_a = n_cart_a+l_A+2 345 n_cartp_b = n_cart_b+l_B+2 346 n_cartm_a = n_cart_a-l_A-1 347 n_cartm_b = n_cart_b-l_B-1 348 n_intpp = n_cartp_a*n_cartp_b*n_cont_ab 349 n_intpm = n_cartm_a*n_cartp_b*n_cont_ab 350 n_intmp = n_cartp_a*n_cartm_b*n_cont_ab 351 n_intmm = n_cartm_a*n_cartm_b*n_cont_ab 352 n_allp_b = n_cartp_b*n_cont_b 353 n_allp_a = n_cartp_a*n_cont_a 354 n_allm_b = n_cartm_b*n_cont_b 355 n_allm_a = n_cartm_a*n_cont_a 356 i_xca = 1 357 i_xcb = i_xca+n_prim_A*n_cont_A 358 i_pp = i_xcb+n_prim_B*n_cont_B 359c i_pm = i_pp+n_intpp*3*nat 360c i_mp = i_pm+n_intpm*3*nat 361c i_mm = i_mp+n_intmp*3*nat 362c i_scr = max(i_xca+n_ab3at*ntyp,i_mm+n_intmm*3*nat) 363 i_pm = i_pp+n_intpp*3*3 364 i_mp = i_pm+n_intpm*3*3 365 i_mm = i_mp+n_intmp*3*3 366 i_scr = max(i_xca+n_ab3at*ntyp,i_mm+n_intmm*3*3) 367 memscr = lscr-i_scr+1 368 369 if (debug_addresses) then 370 write (LuOut,*) 'n_cartp_a',n_cartp_a 371 write (LuOut,*) 'n_cartp_b',n_cartp_b 372 write (LuOut,*) 'n_cartm_a',n_cartm_a 373 write (LuOut,*) 'n_cartm_b',n_cartm_b 374 write (LuOut,*) 'n_intpp',n_intpp 375 write (LuOut,*) 'n_intpm',n_intpm 376 write (LuOut,*) 'n_intmp',n_intmp 377 write (LuOut,*) 'n_intmm',n_intmm 378 write (LuOut,*) 'i_xca,i_xcb',i_xca,i_xcb 379 write (LuOut,*) 'i_pp,i_pm,i_mp,i_mm',i_pp,i_pm,i_mp,i_mm 380 write (LuOut,*) 'i_scr',i_scr 381 call util_flush(LuOut) 382 end if 383* 384* Set up coefficients multiplied by exponents 385* 386 memscr = lscr-i_scr+1 387 if (.not.DryRun) then 388 if (memscr .lt. 0) call errquit ( 389 & 'Insufficient scratch memory in rel_oneld',99, MEM_ERR) 390 k = i_xca-1 391 do j = 1,n_cont_a 392 do i = 1,n_prim_A 393 scr(k+i) = zeta_A(i)*coefS_A(i,j) 394 end do 395 k = k+n_prim_A 396 end do 397 k = i_xcb-1 398 do j = 1,n_cont_B 399 do i = 1,n_prim_B 400 scr(k+i) = zeta_B(i)*coefS_B(i,j) 401 end do 402 k = k+n_prim_B 403 end do 404 end if 405c doS = .false. 406c doT = .false. 407 doV = .true. 408* 409* Calculate integrals for l_A+1, l_B+1 410* 411 if (do_nw) then 412 call hf1d_cosmo( 413 & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1,ictr_A, 414 & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1,ictr_B, 415 & Cxyz,zan,scr(i_pp),n_intpp, 416 & canAB,DryRun,scr(i_scr),memscr) 417 else 418 call errquit("rel_oneld_cosmo: need NW",0,CAPMIS_ERR) 419c call hnd_stvintd( 420c & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1,ictr_A, 421c & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1,ictr_B, 422c & Cxyz,zan,nat,S,T,scr(i_pp),n_intpp,doS,doT,doV, 423c & scr(i_scr),memscr) 424 end if 425 if (debug_arrays) then 426 k = i_pp 427 do i = 1,3 428 do l = 1,3 429 write (LuOut,'(//1x,2A,I3)') xyz(l), 430 & ' component of derivatives for center',i 431 call ecp_matpr(scr(k),1,n_allp_b,1,n_allp_a, 432 & 1,n_allp_b,1,n_allp_a,'l_A+1,l_B+1','E',120,6) 433 k = k+n_ab 434 call util_flush(LuOut) 435 end do 436 end do 437 end if 438 if (DryRun) then 439 max_mem = max(max_mem,i_scr+memscr-1) 440 memscr = lscr-i_scr+1 441 end if 442* 443* Calculate integrals for l_A-1, l_B+1 444* 445 if (l_A .gt. 0) then 446 if (do_nw) then 447 call hf1d_cosmo( 448 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1,ictr_A, 449 & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1,ictr_B, 450 & Cxyz,zan,scr(i_pm),n_intpm, 451 & canAB,DryRun,scr(i_scr),memscr) 452 else 453 call errquit("rel_oneld_cosmo: need NW",0,CAPMIS_ERR) 454c call hnd_stvintd( 455c & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1,ictr_A, 456c & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1,ictr_B, 457c & Cxyz,zan,nat,S,T,scr(i_pm),n_intpm,doS,doT,doV, 458c & scr(i_scr),memscr) 459 end if 460 if (DryRun) then 461 max_mem = max(max_mem,i_scr+memscr-1) 462 memscr = lscr-i_scr+1 463 end if 464 if (debug_arrays) then 465 k = i_pm 466 do i = 1,3 467 do l = 1,3 468 write (LuOut,'(//1x,2A,I3)') xyz(l), 469 & ' component of derivatives for center',i 470 call ecp_matpr(scr(k),1,n_allp_b,1,n_allm_a, 471 & 1,n_allp_b,1,n_allm_a,'l_A-1,l_B+1','E',120,6) 472 k = k+n_ab 473 call util_flush(LuOut) 474 end do 475 end do 476 end if 477 end if 478* 479* Calculate integrals for l_A+1, l_B-1 480* 481 if (l_B .gt. 0) then 482 if (do_nw) then 483 call hf1d_cosmo( 484 & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1,ictr_A, 485 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1,ictr_B, 486 & Cxyz,zan,scr(i_mp),n_intmp, 487 & canAB,DryRun,scr(i_scr),memscr) 488 else 489 call errquit("rel_oneld_cosmo: need NW",0,CAPMIS_ERR) 490c call hnd_stvintd( 491c & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1,ictr_A, 492c & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1,ictr_B, 493c & Cxyz,zan,nat,S,T,scr(i_mp),n_intmp,doS,doT,doV, 494c & scr(i_scr),memscr) 495 end if 496 if (debug_arrays) then 497 k = i_mp 498 do i = 1,3 499 do l = 1,3 500 write (LuOut,'(//1x,2A,I3)') xyz(l), 501 & ' component of derivatives for center',i 502 call ecp_matpr(scr(k),1,n_allm_b,1,n_allp_a, 503 & 1,n_allm_b,1,n_allp_a,'l_A+1,l_B-1','E',120,6) 504 k = k+n_ab 505 call util_flush(LuOut) 506 end do 507 end do 508 end if 509 if (DryRun) then 510 max_mem = max(max_mem,i_scr+memscr-1) 511 memscr = lscr-i_scr+1 512 end if 513* 514* Calculate integrals for l_A-1, l_B-1 515* 516 if (l_A .gt. 0) then 517 if (do_nw) then 518 call hf1d_cosmo( 519 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1,ictr_A, 520 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1,ictr_B, 521 & Cxyz,zan,scr(i_mm),n_intmm, 522 & canAB,DryRun,scr(i_scr),memscr) 523 else 524 call errquit("rel_oneld_cosmo: need NW",0,CAPMIS_ERR) 525c call hnd_stvintd( 526c & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1,ictr_A, 527c & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1,ictr_B, 528c & Cxyz,zan,nat,S,T,scr(i_mm),n_intmm,doS,doT,doV, 529c & scr(i_scr),memscr) 530 end if 531 if (debug_arrays) then 532 k = i_mm 533 do i = 1,3 534 do l = 1,3 535 write (LuOut,'(//1x,2A,I3)') xyz(l), 536 & ' component of derivatives for center',i 537 call ecp_matpr(scr(k),1,n_allm_b,1,n_allm_a, 538 & 1,n_allm_b,1,n_allm_a,'l_A-1,l_B-1','E',120,6) 539 k = k+n_ab 540 call util_flush(LuOut) 541 end do 542 end do 543 end if 544 if (DryRun) then 545 max_mem = max(max_mem,i_scr+memscr-1) 546 memscr = lscr-i_scr+1 547 end if 548 end if 549 end if 550* 551* Compute the relativistic potential energy integrals 552* 553 call rel_pot (scr(i_pp),scr(i_pm),scr(i_mp),scr(i_mm), 554 & scr,n_ab3at,ntyp, 555 & l_A,n_cartp_a,n_cart_a,n_cartm_a,n_cont_A*3*3, 556 & l_B,n_cartp_b,n_cart_b,n_cartm_b,n_cont_B, 557 & DryRun,scr(i_scr),memscr,ibug/10) 558 if (DryRun) then 559 max_mem = max(max_mem,i_scr+memscr-1) 560 lscr = max_mem 561 else 562 i = 1 563 do j = 1,ntyp 564 if (debug_arrays) then 565 write (LuOut,'(//2A)') pot_type(j),' potential' 566 k = i 567 do m = 1,3 568 do l = 1,3 569 write (LuOut,'(//1x,2A,I3)') xyz(l), 570 & ' component of derivatives for center',m 571 call ecp_matpr(scr(k),1,n_all_b,1,n_all_a, 572 & 1,n_all_b,1,n_all_a,'SS potential','E',120,6) 573 k = k+n_ab 574 call util_flush(LuOut) 575 end do 576 end do 577 end if 578 call daxpy (n_ab3at,qalsq,scr(i),1,V(1,j),1) 579 if (debug_arrays) then 580 k = 1 581 do m = 1,3 582 do l = 1,3 583 write (LuOut,'(1x,2A,I3)') xyz(l), 584 & ' component of derivatives for center',m 585 call ecp_matpr(scr(k),1,n_all_b,1,n_all_a,1,n_all_b, 586 & 1,n_all_a,'Relativistic potential','E',120,6) 587 call util_flush(LuOut) 588 k = k+n_ab 589 end do 590 end do 591 end if 592 i = i+n_ab3at 593 end do 594 end if 595 if (debug_gen) write (LuOut,*) 'Exiting rel_oneld' 596* 597 return 598 end 599