1C> \ingroup nwint 2C> @{ 3C> 4C> \brief Compute the electron - nuclear/Bq integral derivatives 5C> 6C> Compute derivatives of the integrals 7C> \f{eqnarray*}{ 8C> \sum_C(\mu|r_C1^{-1}|\rho) = \sum_C\int_{-\infty}^{\infty} 9C> \chi_\mu(R_A;r_1)\chi_\rho(R_B;r_1)|r_1-R_C|^{-1}q_C dr 10C> \f} 11C> These integrals are differentiated with respect to the position of 12C> centers \f$ C \f$. These centers may include both atoms and point charges, 13C> i.e. Bq centers. 14C> 15C> The positions of the centers \f$ C \f$ are picked up from the current 16C> geometry and Bq data structure. 17C> 18C> On output the integral derivatives are stored as if the buffer was allocated 19C> as `H1a(nint,ncoord,ncenter)` where `nint` is the number of integrals in 20C> the shell pair, `ncoord` is the number of Cartesian coordinates, and 21C> `ncenter` is the number of centers. 22C> 23 subroutine intd_1epot(i_basis,ish,j_basis,jsh,lscr,scr, 24 & lH1a,H1a) 25C $Id$ 26 implicit none 27#include "stdio.fh" 28#include "errquit.fh" 29#include "nwc_const.fh" 30#include "basP.fh" 31#include "basdeclsP.fh" 32#include "geomP.fh" 33#include "geobasmapP.fh" 34c 35c layer routine to compute the derivative 1 electron hamiltonian integrals 36c for shells/contractions ish,jsh 37c 38c Order is... nint*3*nat (3=> xyz, nat=number of atoms) 39c 40c / | 41c | nint,d <ij> | 42c | --------------| 43c \ d[idatom(1),x]| 44c | 45c nint,d <ij> | 46c --------------| 47c d[idatom(1),y]| 48c | 49c nint,d <ij> | 50c --------------| 51c d[idatom(1),z]| 52c | 53c nint,d <ij> | 54c --------------| 55c d[idatom(2),x]| 56c | 57c nint,d <ij> | 58c --------------| 59c d[idatom(2),y]| 60c | 61c nint,d <ij> | 62c -------------- | 63c d[idatom(2),z] | 64c 65c . . . 66c | 67c nint,d <ij> | 68c --------------| 69c d[idatom(nat),x]| 70c | 71c nint,d <ij> | 72c --------------| 73c d[idatom(nat),y]| 74c \ 75c nint,d <ij> | 76c -------------- | 77c d[idatom(nat),z]/ 78c 79c::functions 80 integer int_nint_cart 81 external int_nint_cart 82c::passed 83 integer i_basis !< [Input] ish basis set handle 84 integer ish !< [Input] `i` contraction index 85 integer j_basis !< [Input] jsh basis set handle 86 integer jsh !< [Input] `j` contraction index 87 integer lscr !< [Input] length of scratch space 88 integer lH1a !< [Input] number of h1 integral derivatives in shells ish and jsh 89c ! NOTE: nint*3 integral derivatives returned per unique center 90 double precision scr(lscr) !< [Input] scratch array 91 double precision H1a(*) !< [Output] derivative integrals 92c 93c::local 94 integer nint, offset, scrsize, nat 95c 96 nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset))) 97c 98 nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0) 99 if (nint*3*nat.gt.lH1a) then 100 write(luout,*) 'nint*3*nat = ',nint*3*nat 101 write(luout,*) 'lH1a = ',lH1a 102 call errquit('intd_1eh1: nint>lH1a error',911, INT_ERR) 103 endif 104c 105 offset = nint*3*2 ! scratch for Ta array in intd_1eh1P routine 106 scrsize = lscr - offset ! new scratch array size 107 offset = offset + 1 ! increment for passing to intd_1eh1P 108c 109 call intd_1epotP(i_basis,ish,j_basis,jsh, 110cc AJL/Begin/SPIN ECPs 111 & scrsize,scr(offset),nint,H1a,scr,1) 112c ^ ECP Channel (ALPHA) 113c 114 end 115 subroutine intd_1epotP(i_basis,ish,j_basis,jsh,lscr,scr, 116 & nint,H1a,Ta,ecp_channel) 117cc AJL/End 118 implicit none 119#include "stdio.fh" 120#include "errquit.fh" 121#include "apiP.fh" 122#include "nwc_const.fh" 123#include "int_nbf.fh" 124#include "basP.fh" 125#include "basdeclsP.fh" 126#include "geomP.fh" 127#include "geom.fh" 128#include "geobasmapP.fh" 129#include "mafdecls.fh" 130#include "bas_exndcf_dec.fh" 131#include "bas_ibs_dec.fh" 132#include "rel_nwc.fh" 133c::external subroutines used 134c... errquit 135c::functions 136 logical cando_hnd_1e 137 logical cando_nw 138 external cando_hnd_1e 139 external cando_nw 140c::passed 141 integer i_basis !< [Input] ish basis set handle 142 integer ish !< [Input] `i` contraction index 143 integer j_basis !< [Input] jsh basis set handle 144 integer jsh !< [Input] `j` contraction index 145 integer lscr !< [Input] length of scratch space 146 integer nint !< [Input] number of integrals in shells ish and jsh 147c ! NOTE: nint*3 integral derivatives returned per unique center 148 double precision scr(lscr) !< [Input] scratch array 149 double precision H1a(nint,3,*) !< [Output] derivative integrals (nint,3,n_atoms) 150 double precision Ta(nint,3,2) !< [Scratch] space for kinetic integrals 151cc AJL/Begin/SPIN ECPs 152 integer ecp_channel !< [Input] ECP Channel (1 = Alpha, 2 = Beta) 153cc AJL/End 154c::local 155 logical doT 156 integer ucont,uconts 157 integer ibas,iatom,inp,igen,iexp,icf,itype,igeom,isbas,icfS 158 integer jbas,jatom,jnp,jgen,jexp,jcf,jtype,jgeom,jsbas,jcfS 159 integer nat 160 integer nintV 161 integer offset 162c 163 logical any_spherical 164 logical orel, oirel, ojrel, oNR 165 logical ohnd_ok, onw_ok 166 integer i_nbf_x, j_nbf_x 167 integer i_nbf_s, j_nbf_s 168 integer nint_x, nint_s 169 integer zatom, zyx 170 integer lbas, sbas, abas 171c 172 integer ib,bq_ncent,nb 173 integer i_qbq,i_cbq 174 integer i_qbq0,i_cbq0,i_xbq0 175 integer h_qbq0,h_cbq0,h_xbq0 176c 177#include "bas_exndcf_sfn.fh" 178#include "bas_ibs_sfn.fh" 179c 180c check if gencon/sp shells 181c 182 call int_nogencont_check(i_basis,'intd_1eh1P:i_basis') 183 call int_nogencont_check(j_basis,'intd_1eh1P:j_basis') 184 call int_nospshell_check(i_basis,'intd_1eh1P:i_basis') 185 call int_nospshell_check(j_basis,'intd_1eh1P:j_basis') 186c 187 ibas = i_basis + BASIS_HANDLE_OFFSET 188 jbas = j_basis + BASIS_HANDLE_OFFSET 189 oNR = .true. 190 oirel = .false. 191 ojrel = .false. 192 orel = .false. 193c 194 if (dyall_mod_dir) then 195c 196c get basis set handles; relativistic integral option only valid 197c if both ibas and jbas are the ao basis. 198c 199 lbas = lc_bsh + BASIS_HANDLE_OFFSET 200 sbas = sc_bsh + BASIS_HANDLE_OFFSET 201 abas = ao_bsh + BASIS_HANDLE_OFFSET 202 orel = ibas .eq. abas .and. jbas .eq. abas 203 end if 204c 205c i shell 206c 207 ucont = (sf_ibs_cn2ucn(ish,ibas)) 208c 209c check for relativistic shell 210c 211 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then 212 oirel = .true. 213 isbas = sbas 214 uconts = ao_to_ls(ucont) 215 if (uconts .eq. 0) call errquit ( 216 & 'intd_1eh1: no relativistic pointer',911, INT_ERR) 217 if (nesc_1e_approx) then 218 ibas = lbas 219 ucont = uconts 220 end if 221 else 222 uconts = ucont 223 isbas = ibas 224 end if 225c 226 inp = infbs_cont(CONT_NPRIM,ucont,ibas) 227 igen = infbs_cont(CONT_NGEN,ucont,ibas) 228 iexp = infbs_cont(CONT_IEXP,ucont,ibas) 229 icf = infbs_cont(CONT_ICFP,ucont,ibas) 230 itype = infbs_cont(CONT_TYPE,ucont,ibas) 231 igeom = ibs_geom(ibas) 232 iatom = (sf_ibs_cn2ce(ish,ibas)) 233 icfS = infbs_cont(CONT_ICFP ,uconts,isbas) 234c 235c j shell 236c 237 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 238c 239c check for relativistic shell 240c 241 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then 242 ojrel = .true. 243 jsbas = sbas 244 uconts = ao_to_ls(ucont) 245 if (uconts .eq. 0) call errquit ( 246 & 'intd_1eh1: no relativistic pointer',911, INT_ERR) 247 if (nesc_1e_approx) then 248 jbas = lbas 249 ucont = uconts 250 end if 251 else 252 uconts = ucont 253 jsbas = jbas 254 end if 255c 256 jnp = infbs_cont(CONT_NPRIM,ucont,jbas) 257 jgen = infbs_cont(CONT_NGEN,ucont,jbas) 258 jexp = infbs_cont(CONT_IEXP,ucont,jbas) 259 jcf = infbs_cont(CONT_ICFP,ucont,jbas) 260 jtype = infbs_cont(CONT_TYPE,ucont,jbas) 261 jgeom = ibs_geom(jbas) 262 jatom = (sf_ibs_cn2ce(jsh,jbas)) 263 jcfS = infbs_cont(CONT_ICFP ,uconts,jsbas) 264c 265 oNR = .not.(oirel.and.ojrel) 266 orel = oirel.or.ojrel 267c 268 if (igeom.ne.jgeom) then 269 write(luout,*)'intd_1eh1P.F: two different geometries for', 270 & ' derivatives?' 271 call errquit('intd_1eh1P: geom error ',911, INT_ERR) 272 endif 273c 274 if (iatom.eq.jatom) then 275 doT = .false. 276 else 277 doT = .true. 278 endif 279c 280 ohnd_ok = cando_hnd_1e(i_basis,ish,0) 281 & .and. cando_hnd_1e(j_basis,jsh,0) 282 & .and. (.not.geom_any_finuc (igeom)) 283 & .and. (.not.geom_any_finuc (jgeom)) 284 onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0) 285c 286c get external charges here (MV) 287c ------------------------------ 288 if(.not.geom_extbq_on()) 289 > call errquit('int_1eefc:no active bqs',0,0) 290 bq_ncent = geom_extbq_ncenter() 291 i_cbq = geom_extbq_coord() 292 i_qbq = geom_extbq_charge() 293 nb = bq_ncent+ncenter(igeom) 294 295 if(.not.ma_push_get(mt_dbl,3*nb,'cbq',h_cbq0,i_cbq0)) 296 + call errquit( 'intd1_epot',0,0) 297 if(.not.ma_push_get(mt_dbl,nb,'qbq',h_qbq0,i_qbq0)) 298 + call errquit( 'intd1_epot',0,0) 299 if(.not.ma_push_get(mt_dbl,nb,'xbq',h_xbq0,i_xbq0)) 300 + call errquit( 'intd1_epot',0,0) 301 302 call dfill(nb,0.d0,dbl_mb(i_xbq0),1) 303 304 call dcopy(3*ncenter(igeom),coords(1,1,igeom),1, 305 > dbl_mb(i_cbq0),1) 306 call dcopy(3*bq_ncent,dbl_mb(i_cbq),1, 307 > dbl_mb(i_cbq0+3*ncenter(igeom)),1) 308 call dcopy(ncenter(igeom),charge(1,igeom),1, 309 > dbl_mb(i_qbq0),1) 310 call dcopy(bq_ncent,dbl_mb(i_qbq),1, 311 > dbl_mb(i_qbq0+ncenter(igeom)),1) 312 call dcopy(ncenter(igeom),geom_invnucexp(1,igeom),1, 313 > dbl_mb(i_xbq0),1) 314 315c 316 if (orel) then 317 call rel_oneld ( 318 & coords(1,iatom,igeom), 319 & dbl_mb(mb_exndcf(iexp,ibas)), 320 & dbl_mb(mb_exndcf(icf,ibas)), 321 & dbl_mb(mb_exndcf(icfS,isbas)),inp,igen,itype,iatom, 322 & coords(1,jatom,jgeom), 323 & dbl_mb(mb_exndcf(jexp,jbas)), 324 & dbl_mb(mb_exndcf(jcf,jbas)), 325 & dbl_mb(mb_exndcf(jcfS,jsbas)),jnp,jgen,jtype,jatom, 326c & coords(1,1,igeom),charge(1,igeom), 327c & geom_invnucexp(1,igeom),ncenter(igeom), 328 & dbl_mb(i_cbq0),dbl_mb(i_qbq0), 329 & dbl_mb(i_xbq0),nb, 330c........................ doS doT doV canAB 331 & scr,Ta,H1a,nint,.false.,doT,.true.,.false.,onw_ok, 332c........... nonrel dryrun 333 & ohnd_ok,oNR,.false.,scr,lscr,rel_dbg,rel_typ) 334 else if (onw_ok) then 335 call hf1d( 336 & coords(1,iatom,igeom), 337 & dbl_mb(mb_exndcf(iexp,ibas)), 338 & dbl_mb(mb_exndcf(icf,ibas)), 339 & inp,igen,itype,iatom, 340c 341 & coords(1,jatom,jgeom), 342 & dbl_mb(mb_exndcf(jexp,jbas)), 343 & dbl_mb(mb_exndcf(jcf,jbas)), 344 & jnp,jgen,jtype,jatom, 345c 346c & coords(1,1,igeom),charge(1,igeom), 347c & geom_invnucexp(1,igeom),ncenter(igeom), 348 & dbl_mb(i_cbq0),dbl_mb(i_qbq0), 349 & dbl_mb(i_xbq0),nb, 350 & scr,Ta,H1a,nint, 351c..............overlap, k-e, pot-e, canab, dryrun 352 & .false., doT, .true., .false., .false., 353 & scr,lscr) 354 elseif (ohnd_ok) then 355 call hnd_stvintd( 356 & coords(1,iatom,igeom), 357 & dbl_mb(mb_exndcf(iexp,ibas)), 358 & dbl_mb(mb_exndcf(icf,ibas)), 359 & inp,igen,itype,iatom, 360c 361 & coords(1,jatom,jgeom), 362 & dbl_mb(mb_exndcf(jexp,jbas)), 363 & dbl_mb(mb_exndcf(jcf,jbas)), 364 & jnp,jgen,jtype,jatom, 365c 366c & coords(1,1,igeom),charge(1,igeom),ncenter(igeom), 367 & dbl_mb(i_cbq),dbl_mb(i_qbq), 368 & bq_ncent, 369 & scr,Ta,H1a,nint, 370c............overlap, k-e, pot-e, 371 & .false., doT, .true., 372 & scr,lscr) 373 else 374 call errquit('intd_1eh1: could not do hnd or nw integrals', 375 & 0, INT_ERR) 376 endif 377c 378 if(.not.ma_pop_stack(h_xbq0)) 379 + call errquit( 'intd1_epot',0,0) 380 if(.not.ma_pop_stack(h_qbq0)) 381 + call errquit( 'intd1_epot',0,0) 382 if(.not.ma_pop_stack(h_cbq0)) 383 + call errquit( 'intd1_epot',0,0) 384 385c 386c if needed add in Ta derivative integrals 387c 388 if (doT) then 389 call daxpy(nint*3,1.0d00,Ta(1,1,1),1,H1a(1,1,iatom),1) 390 call daxpy(nint*3,1.0d00,Ta(1,1,2),1,H1a(1,1,jatom),1) 391 endif 392c 393c check for ecp 394c 395* 396* this should move to hf1dsp when sp is enabled. 397* 398 nat = ncenter(igeom) ! needed for both ecp and spherical 399 if (any_ecp) then 400 nintV = int_nbf_x(itype)*int_nbf_x(jtype) 401 offset = nintV*3*nat + 1 402* write(luout,*)' lscr to ecp_hf1:',(lscr-nintV) 403* call dcopy(nintV*3*nat,0.0d00,0,scr,1) ! buffer zeroed in ecp_gradient 404c 405cc AJL/Begin/SPIN ECPs 406 if (ecp_channel.eq.1) then 407 call intd_ecp_hf1( 408 & coords(1,iatom,igeom), 409 & dbl_mb(mb_exndcf(iexp,ibas)), 410 & dbl_mb(mb_exndcf(icf,ibas)), 411 & inp,igen,itype,iatom, 412c 413 & coords(1,jatom,jgeom), 414 & dbl_mb(mb_exndcf(jexp,jbas)), 415 & dbl_mb(mb_exndcf(jcf,jbas)), 416 & jnp,jgen,jtype,jatom, 417c 418 & scr,nintV,nat, 419 & scr(offset),(lscr-offset-1), 420 & .false.) 421c 422 else 423 call intd_ecp_hf1_beta( 424 & coords(1,iatom,igeom), 425 & dbl_mb(mb_exndcf(iexp,ibas)), 426 & dbl_mb(mb_exndcf(icf,ibas)), 427 & inp,igen,itype,iatom, 428c 429 & coords(1,jatom,jgeom), 430 & dbl_mb(mb_exndcf(jexp,jbas)), 431 & dbl_mb(mb_exndcf(jcf,jbas)), 432 & jnp,jgen,jtype,jatom, 433c 434 & scr,nintV,nat, 435 & scr(offset),(lscr-offset-1), 436 & .false.) 437c 438 endif 439*... sum ecp into derivative H1 block 440 call daxpy(nintV*3*nat,1.0d00,scr,1,H1a,1) 441cc AJL/End 442 endif 443c 444* H1a now has the cartesian integral block (jlo:jhi,ilo:ihi) 445* 446 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 447 if (.not.any_spherical) return 448c 449c ... reset general contractions for sp shells to 1 since they are handled 450c as a block of 4. 451c 452 if (itype.eq.-1) igen = 1 453 if (jtype.eq.-1) jgen = 1 454c 455 if (bas_spherical(ibas).and.bas_spherical(jbas)) then 456*... transform both i and j integrals 457 i_nbf_x = int_nbf_x(Itype) 458 i_nbf_s = int_nbf_s(Itype) 459 j_nbf_x = int_nbf_x(Jtype) 460 j_nbf_s = int_nbf_s(Jtype) 461c 462 do zatom = 1,nat 463 do zyx = 1,3 464 call spcart_tran1e(H1a(1,zyx,zatom),scr, 465 & j_nbf_x,i_nbf_x,Jtype,jgen, 466 & j_nbf_s,i_nbf_s,Itype,igen, 467 & .false.) 468 enddo 469 enddo 470 else if (bas_spherical(ibas)) then 471*.. transform on i component 472 i_nbf_x = int_nbf_x(Itype) 473 i_nbf_s = int_nbf_s(Itype) 474 j_nbf_x = int_nbf_x(Jtype) 475 j_nbf_s = j_nbf_x 476 do zatom = 1,nat 477 do zyx = 1,3 478 call spcart_tran1e(H1a(1,zyx,zatom),scr, 479 & j_nbf_x,i_nbf_x,0,jgen, 480 & j_nbf_s,i_nbf_s,Itype,igen, 481 & .false.) 482 enddo 483 enddo 484 else if (bas_spherical(jbas)) then 485*.. transform on j component 486 i_nbf_x = int_nbf_x(Itype) 487 i_nbf_s = i_nbf_x 488 j_nbf_x = int_nbf_x(Jtype) 489 j_nbf_s = int_nbf_s(Jtype) 490 do zatom = 1,nat 491 do zyx = 1,3 492 call spcart_tran1e(H1a(1,zyx,zatom),scr, 493 & j_nbf_x,i_nbf_x,Jtype,jgen, 494 & j_nbf_s,i_nbf_s,0,igen, 495 & .false.) 496 enddo 497 enddo 498 else 499 call errquit( 500 & 'int_1eh1: should never reach transform blocked else', 501 & 911, INT_ERR) 502 endif 503c 504c now shuffle transformed buffers to contiguous space 505c 506 nint_x = i_nbf_x*j_nbf_x 507 nint_s = i_nbf_s*j_nbf_s 508 if (nint_s.gt.nint_x) then 509 call errquit 510 & ('intd_1eh1: nint_s >.nint_x diff=',(nint_s-nint_x), 511 & INT_ERR) 512 elseif (nint_s.eq.nint_x) then 513 return 514 else 515 call int_c2s_mv 516 & (H1a,nint_x,nint_s,(3*nat),scr,lscr,'intd_1eh1') 517 endif 518c 519 520 end 521cc AJL/Begin/SPIN ECPs 522cc This is a wrapper, identical to intd_1epot except we have separated out 523cc For alpha and beta channels 524 subroutine intd_1epot_beta(i_basis,ish,j_basis,jsh,lscr,scr, 525 & lH1a,H1a) 526 implicit none 527#include "stdio.fh" 528#include "errquit.fh" 529#include "nwc_const.fh" 530#include "basP.fh" 531#include "basdeclsP.fh" 532#include "geomP.fh" 533#include "geobasmapP.fh" 534 integer int_nint_cart 535 external int_nint_cart 536c::passed 537 integer i_basis !< [Input] ish basis set handle 538 integer ish !< [Input] `i` contraction index 539 integer j_basis !< [Input] jsh basis set handle 540 integer jsh !< [Input] `j` contraction index 541 integer lscr !< [Input] length of scratch space 542 integer lH1a !< [Input] number of h1 integral derivatives in shells ish and jsh 543c ! NOTE: nint*3 integral derivatives returned per unique center 544 double precision scr(lscr) !< [Input] scratch array 545 double precision H1a(*) !< [Output] derivative integrals 546c::local 547 integer nint, offset, scrsize, nat 548c 549 nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset))) 550c 551 nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0) 552 if (nint*3*nat.gt.lH1a) then 553 write(luout,*) 'nint*3*nat = ',nint*3*nat 554 write(luout,*) 'lH1a = ',lH1a 555 call errquit('intd_1eh1: nint>lH1a error',911, INT_ERR) 556 endif 557c 558 offset = nint*3*2 ! scratch for Ta array in intd_1eh1P routine 559 scrsize = lscr - offset ! new scratch array size 560 offset = offset + 1 ! increment for passing to intd_1eh1P 561c 562 call intd_1epotP(i_basis,ish,j_basis,jsh, 563 & scrsize,scr(offset),nint,H1a,scr,2) 564c ^ ECP Channel (BETA) 565 end 566C> @} 567