1C> \ingroup nwint 2C> @{ 3C> 4C> \brief Compute the integral 2nd derivatives of the 1-electron Hamiltonian 5C> 6C> Compute the integral 2nd derivatives of the 1-electron Hamiltonian 7C> consisting of the kinetic energy integrals and the nuclear attraction 8C> integrals. The kinetic energy integral 2nd derivatives are given by 9C> \f{eqnarray*}{ 10C> -\frac{1}{2}\frac{\partial^2 (\mu|\nabla_1^2|\nu)}{\partial R_x\partial R_y} &=& 11C> -\frac{1}{2}\int \frac{\partial^2 [g_\mu(X_\mu,r_1)\nabla_1^2g_\nu(X_\nu,r_1)]} 12C> {\partial X_x\partial X_y}dr_1 13C> \f} 14C> The nuclear attraction integral 2nd derivatives are given by 15C> \f{eqnarray*}{ 16C> \sum_A\frac{\partial^2 (\mu|Z_AR_A^{-1}|\nu)}{\partial R_x\partial R_y} &=& 17C> \sum_A\int \frac{\partial^2 [g_\mu(X_\mu,r_1)Z_AR^{-1}_{A1}g_\nu(X_\nu,r_1)]} 18C> {\partial X_x\partial R_y}dr_1 19C> \f} 20C> The integral 2nd derivatives are returned in `H1a` in an order that is 21C> equivalent to the declaration `H1a(nint,ncoordu,ncoordv,natom,ncoorda)`, 22C> where `nint` refers to the number of integrals associated with the shell 23C> pair, `ncoordu` refers the number of Cartesian coordinates of the atom 24C> associated shell `ish` likewise `ncoordv` refers to atomic coordinates of the 25C> `jsh` atom, `natom` is the number of atoms and `ncoorda` is the number 26C> of coordinates of each nucleus. 27C> 28 subroutine intdd_1eh1(i_basis,ish,j_basis,jsh,lscr,scr, 29 & lH1a,H1a) 30C $Id$ 31 implicit none 32#include "stdio.fh" 33#include "errquit.fh" 34#include "nwc_const.fh" 35#include "basP.fh" 36#include "basdeclsP.fh" 37#include "geomP.fh" 38#include "geobasmapP.fh" 39c 40c layer routine to compute the derivative 1 electron hamiltonian integrals 41c for shells/contractions ish,jsh 42c 43c Order is... nint*3*nat (3=> xyz, nat=number of atoms) 44c 45c / | 46c | nint,d <ij> | 47c | --------------| 48c \ d[idatom(1),x]| 49c | 50c nint,d <ij> | 51c --------------| 52c d[idatom(1),y]| 53c | 54c nint,d <ij> | 55c --------------| 56c d[idatom(1),z]| 57c | 58c nint,d <ij> | 59c --------------| 60c d[idatom(2),x]| 61c | 62c nint,d <ij> | 63c --------------| 64c d[idatom(2),y]| 65c | 66c nint,d <ij> | 67c -------------- | 68c d[idatom(2),z] | 69c 70c . . . 71c | 72c nint,d <ij> | 73c --------------| 74c d[idatom(nat),x]| 75c | 76c nint,d <ij> | 77c --------------| 78c d[idatom(nat),y]| 79c \ 80c nint,d <ij> | 81c -------------- | 82c d[idatom(nat),z]/ 83c 84c::functions 85 integer int_nint_cart 86 external int_nint_cart 87c::passed 88 integer i_basis !< [Input] ish basis set handle 89 integer ish !< [Input] `i` contraction index 90 integer j_basis !< [Input] jsh basis set handle 91 integer jsh !< [Input] `j` contraction index 92 integer lscr !< [Input] length of scratch space 93 integer lH1a !< [Input] number of h1 integral derivatives in shells ish and jsh 94c ! NOTE: nint*3 integral derivatives returned per unique center 95 double precision scr(lscr) !< [Input] scratch array 96 double precision H1a(*) !< [Output] derivative integrals 97c 98c::local 99 integer nint, offset, nat 100c 101 nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset))) 102c 103 nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0) 104 if (nint*3*3*(nat*3+3).gt.lH1a) then 105 write(luout,*) 'nint*3*3*(nat*3+3) = ',nint*3*3*(nat*3+3) 106 write(luout,*) 'lH1a = ',lH1a 107 call errquit('intdd_1eh1: nint>lH1a error',911, INT_ERR) 108 endif 109c 110 offset = nint*3*3*nat*3 + 1 ! T is held seperately in H1a 111c 112 call intdd_1eh1P(i_basis,ish,j_basis,jsh, 113 & lscr,scr,nint,H1a,H1a(offset),nat) 114c 115 end 116 subroutine intdd_1eh1P(i_basis,ish,j_basis,jsh,lscr,scr, 117 & nint,H1a,Ta,nat) 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 "geobasmapP.fh" 128#include "mafdecls.fh" 129#include "bas_exndcf_dec.fh" 130#include "bas_ibs_dec.fh" 131c::external subroutines used 132c... errquit 133c::functions 134 logical cando_hnd_1edd 135 logical cando_nw 136 external cando_hnd_1edd 137 external cando_nw 138c::passed 139 integer i_basis ! [input] ish basis set handle 140 integer ish ! [input] ``i'' contraction index 141 integer j_basis ! [input] jsh basis set handle 142 integer jsh ! [input] ``j'' contraction index 143 integer lscr ! [input] length of scratch space 144 integer nat ! [input] number of atoms 145 integer nint ! [input] number of integrals in shells ish and jsh 146c ! NOTE: nint*3*3 integral derivatives returned per unique center 147 double precision scr(lscr) ! [input] scratch array 148 double precision H1a(nint,3,3,*) ! [output] derivative integrals (nint,3,3,n_atoms,3) 149 double precision Ta(nint,3,3,3) ! [scratch] space for kinetic integrals 150c::local 151 logical doT 152 integer ucont 153 integer ibas,iatom,inp,igen,iexp,icf,itype,igeom 154 integer jbas,jatom,jnp,jgen,jexp,jcf,jtype,jgeom 155 integer i_nbf_x, j_nbf_x 156 integer i_nbf_s, j_nbf_s 157 integer nint_x, nint_s 158 integer zatom, zyx1, zyx2 159c 160 logical any_spherical 161c 162c Temporary variable that needs to be taken out after testing! 163c 164c integer itemp,jtemp,ktemp,ltemp 165c 166#include "bas_exndcf_sfn.fh" 167#include "bas_ibs_sfn.fh" 168c 169c check if gencon/sp shells 170c 171 call int_nogencont_check(i_basis,'intd_1eh1P:i_basis') 172 call int_nogencont_check(j_basis,'intd_1eh1P:j_basis') 173 call int_nospshell_check(i_basis,'intd_1eh1P:i_basis') 174 call int_nospshell_check(j_basis,'intd_1eh1P:j_basis') 175c 176 ibas = i_basis + BASIS_HANDLE_OFFSET 177 jbas = j_basis + BASIS_HANDLE_OFFSET 178c 179 ucont = (sf_ibs_cn2ucn(ish,ibas)) 180 inp = infbs_cont(CONT_NPRIM,ucont,ibas) 181 igen = infbs_cont(CONT_NGEN,ucont,ibas) 182 iexp = infbs_cont(CONT_IEXP,ucont,ibas) 183 icf = infbs_cont(CONT_ICFP,ucont,ibas) 184 itype = infbs_cont(CONT_TYPE,ucont,ibas) 185 igeom = ibs_geom(ibas) 186 iatom = (sf_ibs_cn2ce(ish,ibas)) 187c 188 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 189 jnp = infbs_cont(CONT_NPRIM,ucont,jbas) 190 jgen = infbs_cont(CONT_NGEN,ucont,jbas) 191 jexp = infbs_cont(CONT_IEXP,ucont,jbas) 192 jcf = infbs_cont(CONT_ICFP,ucont,jbas) 193 jtype = infbs_cont(CONT_TYPE,ucont,jbas) 194 jgeom = ibs_geom(jbas) 195 jatom = (sf_ibs_cn2ce(jsh,jbas)) 196c 197 if (igeom.ne.jgeom) then 198 write(luout,*)'intdd_1eh1P.F: two different geometries for', 199 & ' derivatives?' 200 call errquit('intdd_1eh1P: geom error ',911, GEOM_ERR) 201 endif 202c 203 if (iatom.eq.jatom) then 204 doT = .false. 205 else 206 doT = .true. 207 endif 208c 209 if (cando_hnd_1edd(i_basis,ish,0).and. 210 & cando_hnd_1edd(j_basis,jsh,0)) then 211 call hnd_stvintdd( 212 & coords(1,iatom,igeom), 213 & dbl_mb(mb_exndcf(iexp,ibas)), 214 & dbl_mb(mb_exndcf(icf,ibas)), 215 & inp,igen,itype,iatom, 216c 217 & coords(1,jatom,jgeom), 218 & dbl_mb(mb_exndcf(jexp,jbas)), 219 & dbl_mb(mb_exndcf(jcf,jbas)), 220 & jnp,jgen,jtype,jatom, 221c 222 & coords(1,1,igeom),charge(1,igeom),ncenter(igeom), 223 & scr,Ta,H1a,nint, 224c............overlap, k-e, pot-e, 225 & .false., doT, .true., 226 & scr,lscr) 227 else 228 call errquit('intdd_1eh1: could not do hnd integrals', 229 & 0, INT_ERR) 230 endif 231c 232 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 233 if (.not.any_spherical) return 234c 235c ... reset general contractions for sp shells to 1 since they are handled 236c as a block of 4. 237c 238 if (itype.eq.-1) igen = 1 239 if (jtype.eq.-1) jgen = 1 240c 241 if (bas_spherical(ibas).and.bas_spherical(jbas)) then 242*... transform both i and j integrals 243 i_nbf_x = int_nbf_x(Itype) 244 i_nbf_s = int_nbf_s(Itype) 245 j_nbf_x = int_nbf_x(Jtype) 246 j_nbf_s = int_nbf_s(Jtype) 247c 248 do zatom = 1,nat*3 249 do zyx2 = 1,3 250 do zyx1 = 1,3 251 call spcart_tran1e(H1a(1,zyx1,zyx2,zatom),scr, 252 & j_nbf_x,i_nbf_x,Jtype,jgen, 253 & j_nbf_s,i_nbf_s,Itype,igen, 254 & .false.) 255 enddo 256 enddo 257 enddo 258 do zatom = 1,3 259 do zyx2 = 1,3 260 do zyx1 = 1,3 261 call spcart_tran1e(Ta(1,zyx1,zyx2,zatom),scr, 262 & j_nbf_x,i_nbf_x,Jtype,jgen, 263 & j_nbf_s,i_nbf_s,Itype,igen, 264 & .false.) 265 enddo 266 enddo 267 enddo 268 else if (bas_spherical(ibas)) then 269*.. transform on i component 270 i_nbf_x = int_nbf_x(Itype) 271 i_nbf_s = int_nbf_s(Itype) 272 j_nbf_x = int_nbf_x(Jtype) 273 j_nbf_s = j_nbf_x 274 do zatom = 1,nat*3 275 do zyx2 = 1,3 276 do zyx1 = 1,3 277 call spcart_tran1e(H1a(1,zyx1,zyx2,zatom),scr, 278 & j_nbf_x,i_nbf_x,0,jgen, 279 & j_nbf_s,i_nbf_s,Itype,igen, 280 & .false.) 281 enddo 282 enddo 283 enddo 284 do zatom = 1,3 285 do zyx2 = 1,3 286 do zyx1 = 1,3 287 call spcart_tran1e(Ta(1,zyx1,zyx2,zatom),scr, 288 & j_nbf_x,i_nbf_x,0,jgen, 289 & j_nbf_s,i_nbf_s,Itype,igen, 290 & .false.) 291 enddo 292 enddo 293 enddo 294 else if (bas_spherical(jbas)) then 295*.. transform on j component 296 i_nbf_x = int_nbf_x(Itype) 297 i_nbf_s = i_nbf_x 298 j_nbf_x = int_nbf_x(Jtype) 299 j_nbf_s = int_nbf_s(Jtype) 300 do zatom = 1,nat*3 301 do zyx2 = 1,3 302 do zyx1 = 1,3 303 call spcart_tran1e(H1a(1,zyx1,zyx2,zatom),scr, 304 & j_nbf_x,i_nbf_x,Jtype,jgen, 305 & j_nbf_s,i_nbf_s,0,igen, 306 & .false.) 307 enddo 308 enddo 309 enddo 310 do zatom = 1,3 311 do zyx2 = 1,3 312 do zyx1 = 1,3 313 call spcart_tran1e(Ta(1,zyx1,zyx2,zatom),scr, 314 & j_nbf_x,i_nbf_x,Jtype,jgen, 315 & j_nbf_s,i_nbf_s,0,igen, 316 & .false.) 317 enddo 318 enddo 319 enddo 320 else 321 call errquit( 322 & 'intdd_1eh1P: cant do sphericals', 323 & 911, INT_ERR) 324 endif 325c 326c now shuffle transformed buffers to contiguous space 327c 328 nint_x = i_nbf_x*j_nbf_x 329 nint_s = i_nbf_s*j_nbf_s 330 if (nint_s.gt.nint_x) then 331 call errquit 332 & ('intdd_1eh1: nint_s >.nint_x diff=',(nint_s-nint_x), 333 & INT_ERR) 334 elseif (nint_s.eq.nint_x) then 335 return 336 else 337 call int_c2s_mv ! do both H1a and Ta at the same time 338 & (H1a,nint_x,nint_s,(27*nat+27),scr,lscr,'intdd_1eh1') 339 endif 340c 341 end 342C> @} 343