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