1c $Id$ 2* 3C> \ingroup nwint 4C> @{ 5C> 6C> \brief Calculate the 4-center 2-electron integrals between GIAO basis functions 7C> 8C> The GIAO basis functions are important in the treatment of magnetic 9C> properties. These properties are not automatically gauge invariant but 10C> by multiplying the usual Atomic Orbitals (AO) with a magnetic field dependent 11C> factor Gauge Invariant Atomic Orbitals (GIAO) can be defined. Using these 12C> basis functions the properties can be cast in an gauge invariant form (see 13C> [1,2]). 14C> 15C> In the presence of a magnetic field a vector potential is generated given by 16C> \f{eqnarray*}{ 17C> \vec{A}_A = \vec{B} \otimes \vec{R}_A 18C> \f} 19C> where \f$ \vec{B} \f$ is the magnetic field, \f$ \vec{R}_A \f$ is the 20C> nuclear position relative to the origin of the systems, 21C> \f$ \vec{A} \f$ the generated vector potential. 22C> A GIAO is now given by 23C> \f{eqnarray*}{ 24C> g^{GIAO}_\mu(X_\mu,r) = g_\mu(X_\mu,r)e^{-i \vec{A}_A\cdot\vec{r}} 25C> \f} 26C> The integrals computed here can now be represented as 27C> \f{eqnarray*}{ 28C> ({\mu}{\rho}|{\nu}{\lambda}) = \int_{-\infty}^{\infty} 29C> g^{GIAO}_{\mu}(X_{\mu},r_{1})g^{GIAO}_{\rho}(X_{\rho},r_{1})\frac{1}{r_{12}} 30C> g^{GIAO}_{\nu}(X_{\nu},r_{2})g^{GIAO}_{\lambda}(X_{\lambda},r_{2})dr_{1}dr_{2} 31C> \f} 32C> 33C> [1] M. Dupuis, 34C> "New integral transforms for molecular properties and application 35C> to a massively parallel GIAO-SCF implementation", 36C> Comp. Phys. Comm. <b>134</b>, 150-166 (2001), DOI: 37C> <a href="https://doi.org/10.1016/S0010-4655(00)00195-8"> 38C> 10.1016/S0010-4655(00)00195-8</a>. 39C> 40C> [2] K. Wolinski, J.F. Hinton, P. Pulay, 41C> "Efficient implementation of the gauge-independent atomic orbital 42C> method for NMR chemical shift calculations", 43C> J. Am. Chem. Soc. <b>112</b>, 8251-8260 (1990), DOI: 44C> <a href="https://doi.org/10.1021/ja00179a005">10.1021/ja00179a005</a>. 45C> 46 subroutine int_giao_2e(brain, ish, jsh, ketin, ksh, lsh, 47 & lscr, scr, leri, eri) 48 implicit none 49c 50c basic api routine to generate a block of two electron integrals 51c eri = <bra_g(ish).bra_g(jsh) | ket_g(ksh).ket_g(lsh)> 52c 53#include "bas.fh" 54#include "errquit.fh" 55#include "nwc_const.fh" 56#include "apiP.fh" 57#include "basP.fh" 58#include "basdeclsP.fh" 59#include "geomP.fh" 60#include "geobasmapP.fh" 61#include "mafdecls.fh" 62#include "bas_exndcf_dec.fh" 63#include "bas_ibs_dec.fh" 64#include "int_nbf.fh" 65#include "stdio.fh" 66#include "rel_nwc.fh" 67#include "util.fh" 68#include "global.fh" 69 common/testdata/timing(20),irepeat 70 double precision timing 71 integer irepeat 72c 73c::external subroutines used 74c errquit 75c::functions 76 logical cando_nw 77 logical cando_sp 78 logical cando_txs 79 external cando_nw 80 external cando_sp 81 external cando_txs 82c:: passed 83 integer brain !< [Input] bra basis set handle 84 integer ish !< [Input] shell/contraction index 85 integer jsh !< [Input] shell/contraction index 86 integer ketin !< [Input] ket basis set handle 87 integer ksh !< [Input] shell/contraction index 88 integer lsh !< [Input] shell/contraction index 89 integer lscr !< [Input] length of scratch array 90 double precision scr(lscr) !< [Scratch] array 91 integer leri !< [In|Output] length of integral array 92 double precision eri(*) !< [Output] 2e4c integrals 93c:: local 94 integer bra, ket 95 integer ab_geom, cd_geom, ucont, ityp 96 integer La, a_prim, a_gen, a_iexp, a_icfp, a_cent, a_icfps 97 integer Lb, b_prim, b_gen, b_iexp, b_icfp, b_cent, b_icfps 98 integer Lc, c_prim, c_gen, c_iexp, c_icfp, c_cent, c_icfps 99 integer Ld, d_prim, d_gen, d_iexp, d_icfp, d_cent, d_icfps 100c 101 double precision roff(3) 102 double precision q4 103 integer nint 104 logical dum_log, do_bra, do_ket 105 logical status_sp, status_nw, status_txs, status_gen 106 integer texas_ang_limit 107 logical ieqj, keql 108 integer sbas, abas, bras, kets 109c 110 logical any_spherical 111 integer a_nbf, b_nbf, c_nbf, d_nbf 112 integer a_nbf_s, b_nbf_s, c_nbf_s, d_nbf_s 113 integer ab_gen, ab_cmp, cd_gen, cd_cmp,i,j 114c 115#include "bas_exndcf_sfn.fh" 116#include "bas_ibs_sfn.fh" 117c 118c timing(10)=timing(10)-util_wallsec() 119 bra = brain + BASIS_HANDLE_OFFSET 120 ket = ketin + BASIS_HANDLE_OFFSET 121 bras = bra 122 kets = ket 123 ab_geom = ibs_geom(bra) 124 cd_geom = ibs_geom(ket) 125 a_cent = (sf_ibs_cn2ce(ish,bra)) 126 b_cent = (sf_ibs_cn2ce(jsh,bra)) 127 c_cent = (sf_ibs_cn2ce(ksh,ket)) 128 d_cent = (sf_ibs_cn2ce(lsh,ket)) 129c 130c 131 any_spherical = bas_spherical(bra).or.bas_spherical(ket) 132c 133 status_sp = cando_sp(brain,ish,jsh).and.cando_sp(ketin,ksh,lsh) 134 if (.not.status_sp) then 135c 136 ieqj = ish.eq.jsh 137 keql = ksh.eq.lsh 138c 139 ucont = sf_ibs_cn2ucn(ish,bra) 140 La = infbs_cont(CONT_TYPE ,ucont,bra) 141 a_prim = infbs_cont(CONT_NPRIM,ucont,bra) 142 a_gen = infbs_cont(CONT_NGEN ,ucont,bra) 143 a_iexp = infbs_cont(CONT_IEXP ,ucont,bra) 144 a_icfp = infbs_cont(CONT_ICFP ,ucont,bra) 145 a_icfps = infbs_cont(CONT_ICFP ,ucont,bras) 146c 147 ucont = sf_ibs_cn2ucn(jsh,bra) 148 Lb = infbs_cont(CONT_TYPE ,ucont,bra) 149 b_prim = infbs_cont(CONT_NPRIM,ucont,bra) 150 b_gen = infbs_cont(CONT_NGEN ,ucont,bra) 151 b_iexp = infbs_cont(CONT_IEXP ,ucont,bra) 152 b_icfp = infbs_cont(CONT_ICFP ,ucont,bra) 153 b_icfps = infbs_cont(CONT_ICFP ,ucont,bras) 154c 155 ucont = sf_ibs_cn2ucn(ksh,ket) 156 Lc = infbs_cont(CONT_TYPE ,ucont,ket) 157 c_prim = infbs_cont(CONT_NPRIM,ucont,ket) 158 c_gen = infbs_cont(CONT_NGEN ,ucont,ket) 159 c_iexp = infbs_cont(CONT_IEXP ,ucont,ket) 160 c_icfp = infbs_cont(CONT_ICFP ,ucont,ket) 161 c_icfps = infbs_cont(CONT_ICFP ,ucont,kets) 162c 163 ucont = sf_ibs_cn2ucn(lsh,ket) 164 Ld = infbs_cont(CONT_TYPE ,ucont,ket) 165 d_prim = infbs_cont(CONT_NPRIM,ucont,ket) 166 d_gen = infbs_cont(CONT_NGEN ,ucont,ket) 167 d_iexp = infbs_cont(CONT_IEXP ,ucont,ket) 168 d_icfp = infbs_cont(CONT_ICFP ,ucont,ket) 169 d_icfps = infbs_cont(CONT_ICFP ,ucont,kets) 170c 171 a_nbf = int_nbf_x(La) 172 b_nbf = int_nbf_x(Lb) 173 c_nbf = int_nbf_x(Lc) 174 d_nbf = int_nbf_x(Ld) 175 leri=a_nbf*b_nbf*c_nbf*d_nbf*a_gen*b_gen*c_gen*d_gen 176 call hnd_giahnd( 177 & coords(1,a_cent,ab_geom), dbl_mb(mb_exndcf(a_iexp,bra)), 178 & dbl_mb(mb_exndcf(a_icfp,bra)), a_prim, a_gen, La, 179 & coords(1,b_cent,ab_geom), dbl_mb(mb_exndcf(b_iexp,bra)), 180 & dbl_mb(mb_exndcf(b_icfp,bra)), b_prim, b_gen, Lb, 181 & coords(1,c_cent,cd_geom), dbl_mb(mb_exndcf(c_iexp,ket)), 182 & dbl_mb(mb_exndcf(c_icfp,ket)), c_prim, c_gen, Lc, 183 & coords(1,d_cent,cd_geom), dbl_mb(mb_exndcf(d_iexp,ket)), 184 & dbl_mb(mb_exndcf(d_icfp,ket)), d_prim,d_gen,Ld, 185 & ieqj, keql, eri, leri, scr, lscr) 186c 187c eri has cartesian block of integrals (llo:lhi,klo:khi,jlo:jhi,ilo:ihi,6) 188c 189 if (any_spherical) then 190 a_nbf_s = int_nbf_s(La) 191 b_nbf_s = int_nbf_s(Lb) 192 c_nbf_s = int_nbf_s(Lc) 193 d_nbf_s = int_nbf_s(Ld) 194 cd_gen = c_gen*d_gen 195 ab_gen = a_gen*b_gen 196 ab_cmp = a_nbf*b_nbf 197 cd_cmp = c_nbf*d_nbf 198 do_bra=bas_spherical(bra) 199 do_ket=bas_spherical(ket) 200 call giao_to_sph(eri,leri,scr,La,Lb,Lc,Ld,a_nbf,b_nbf,c_nbf, 201 & d_nbf,a_nbf_s,b_nbf_s,c_nbf_s,d_nbf_s, 202 & a_gen,b_gen,c_gen,d_gen,ab_gen,ab_cmp, 203 & cd_gen,cd_cmp,do_bra,do_ket) 204 endif 205 else 206 write(luout,*)'int_giao_2e: cannot do sp integrals' 207 write(luout,*)' brain :',brain 208 write(luout,*)' ketin :',ketin 209 write(luout,*)' ish :',ish 210 write(luout,*)' jsh :',jsh 211 write(luout,*)' ksh :',ksh 212 write(luout,*)' lsh :',lsh 213 call errquit('int_giao_2e: fatal error ',0, INT_ERR) 214 endif 215 end 216c 217 subroutine giao_to_sph(eri,leri,scr,La,Lb,Lc,Ld,a_nbf,b_nbf,c_nbf, 218 & d_nbf,a_nbf_s,b_nbf_s,c_nbf_s,d_nbf_s, 219 & a_gen,b_gen,c_gen,d_gen,ab_gen,ab_cmp, 220 & cd_gen,cd_cmp,do_bra,do_ket) 221 implicit none 222c 223 integer leri 224 double precision eri(leri,6), scr(*) 225 integer La,Lb,Lc,Ld,a_nbf,b_nbf,c_nbf,d_nbf 226 integer a_nbf_s,b_nbf_s,c_nbf_s,d_nbf_s,a_gen,b_gen,c_gen,d_gen 227 integer ab_gen,ab_cmp,cd_gen,cd_cmp 228 logical do_bra,do_ket 229c 230 integer ityp 231c 232 ab_cmp = a_nbf_s*b_nbf_s 233 do ityp = 1, 6 234 call spcart_bra2etran(eri(1,ityp),scr, 235 & b_nbf,a_nbf,b_nbf_s,a_nbf_s, 236 & Lb, La, b_gen, a_gen, 237 & cd_gen*cd_cmp,.false.) 238 call spcart_ket2etran(eri(1,ityp),scr, 239 & d_nbf,c_nbf,d_nbf_s,c_nbf_s, 240 & Ld, Lc, d_gen, c_gen, 241 & ab_gen*ab_cmp,.false.) 242 enddo 243 a_nbf = a_nbf_s 244 b_nbf = b_nbf_s 245 c_nbf = c_nbf_s 246 d_nbf = d_nbf_s 247 end 248c 249C> @} 250