1c $Id$ 2* 3C> \ingroup nwint 4C> @{ 5C> \brief Compute any 1-electron integrals 6C> 7C> This is an internal routine that most of the external 1-electron 8C> routines call. This is the actual workhorse routine. 9C> This routine computes the 1 electron integrals \f$S\f$, \f$T\f$, and \f$V\f$: 10C> \f{eqnarray*}{ 11C> S & = & ({\mu}|{\nu}) \\\\ 12C> & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\\\ 13C> T & = & ({\mu}|-\frac{1}{2}\nabla^{2}|{\nu}) \\\\ 14C> & = & -\frac{1}{2}\int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\nabla^{2}(r_{1})\\ 15C> g_{\nu}(X_{\nu},r_{1})dr_{1} \\\\ 16C> V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\\\ 17C> & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 18C> {-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 19C> \f} 20C> 21C> If an ECP is defined then the ECP integral contributions are summed 22C> directly into the V integrals. 23C> 24C> If a relativistic basis is defined then the 1-electron integrals 25C> for the case where both shells are relativistic are modified to 26C> \f{eqnarray*}{ 27C> S & = & ({\mu^L}|{\nu^L}) 28C> - ({\mu^S}|\frac{\alpha^2}{4}{\nabla^{2}}|{\nu^S}) \\\\ 29C> T & = & -\frac{1}{2} ({\mu^L}|{\nabla^{2}}|{\nu^S}) 30C> - \frac{1}{2} ({\mu^S}|{\nabla^{2}}|{\nu^L}) 31C> + \frac{1}{2} ({\mu^S}|{\nabla^{2}}|{\nu^S}) \\\\ 32C> V & = & ({\mu^L}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu^L}) 33C> - \frac{\alpha^2}{4} ({\mu^S}|\nabla\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}\cdot\nabla|{\nu^S}) \\\\ 34C> \f} 35C> 36c:tex-% this is part of the API Standard Integral routines. 37c:tex-\subsection{int\_1estv} 38c:tex-This is an internal routine that most of the external 1 electron 39c:tex-routines call. This is the actual workhorse routine. 40c:tex-This routine computes the 1 electron integrals S, T, and V: 41c:tex-\begin{eqnarray*} 42c:tex-S & = & ({\mu}|{\nu}) \\ 43c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\ 44c:tex-T & = & ({\mu}|-\frac{1}{2}\nabla^{2}|{\nu}) \\ 45c:tex- & = & -\frac{1}{2}\int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\nabla^{2}(r_{1}) 46c:tex- g_{\nu}(X_{\nu},r_{1})dr_{1} \\ 47c:tex-V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\ 48c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 49c:tex-{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} \\ 50c:tex-\end{eqnarray*} 51c:tex- 52c:tex-If an ECP is defined then the ECP integral contributions are summed 53c:tex-directly into the V integrals. 54c:tex- 55c:tex-If a relativistic basis is defined then the one-electron integrals 56c:tex-for the case where both shells are relativistic are modified to 57c:tex-\begin{eqnarray*} 58c:tex-S & = & ({\mu^L}|{\nu^L}) 59c:tex- - ({\mu^S}|\frac{\alpha^2}{4}{\nabla^{2}}|{\nu^S}) \\ 60c:tex-T & = & -\frac{1}{2} ({\mu^L}|{\nabla^{2}}|{\nu^S}) 61c:tex- - \frac{1}{2} ({\mu^S}|{\nabla^{2}}|{\nu^L}) 62c:tex- + \frac{1}{2} ({\mu^S}|{\nabla^{2}}|{\nu^S}) \\ 63c:tex-V & = & ({\mu^L}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu^L}) 64c:tex- - \frac{\alpha^2}{4} ({\mu^S}|\nabla\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}\cdot\nabla|{\nu^S}) \\ 65c:tex-\end{eqnarray*} 66c:tex- 67c:tex- 68c:tex-{\it Syntax:} 69c:tex-\begin{verbatim} 70 subroutine int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V, 71cc AJL/Begin 72 & doS,doT,doV,msg) 73cc AJL/End 74c:tex-\end{verbatim} 75 implicit none 76#include "nwc_const.fh" 77#include "errquit.fh" 78#include "basP.fh" 79#include "basdeclsP.fh" 80#include "geomP.fh" 81#include "geom.fh" 82#include "geobasmapP.fh" 83#include "mafdecls.fh" 84#include "bas_exndcf_dec.fh" 85#include "bas_ibs_dec.fh" 86#include "int_nbf.fh" 87#include "stdio.fh" 88#include "apiP.fh" 89#include "rel_nwc.fh" 90#include "util.fh" 91c::external subroutines used 92c... errquit 93c::functions 94 logical cando_hnd_1e 95 logical cando_nw_1e 96 logical cando_nw 97 logical int_chk_init 98 logical int_chk_sh 99 external int_chk_init 100 external int_chk_sh 101 external cando_hnd_1e 102 external cando_nw_1e 103 external cando_nw 104c::passed 105c:tex-\begin{verbatim} 106 integer i_basis !< [Input] basis set handle for ish 107 integer ish !< [Input] i shell/contraction 108 integer j_basis !< [Input] basis set handle for jsh 109 integer jsh !< [Input] j shell/contraction 110 integer lscr !< [Input] length of scratch array 111 integer lstv !< [Input] length of integral buffer 112 double precision scr(lscr) !< [Scratch] scratch array 113 double precision S(lstv) !< [Output] overlap integrals 114 double precision T(lstv) !< [Output] kinetic energy integrals 115 double precision V(lstv) !< [Output] potential energy integrals 116 logical doS !< [Input] flag for overlap integrals 117 logical doT !< [Input] flag for kinetic energy integrals 118 logical doV !< [Input] flag for potential energy integrals 119c:tex-\end{verbatim} 120c::local 121 logical ohnd_ok, onw_ok, onw1e_ok 122 logical shells_ok, orel, oirel, ojrel, oNR, canAB 123 integer i_geom, j_geom, ibas, jbas, ucont, uconts 124 integer lbas, sbas, abas, isbas, jsbas 125 integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_icfpS 126 integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_icfpS 127c 128 logical any_spherical 129 integer i_nbf_x, j_nbf_x 130 integer i_nbf_s, j_nbf_s 131c 132 integer WarnP 133 save WarnP 134 data WarnP /0/ 135cc AJL/Begin 136cc Soft-coded the value of msg 137 character*(*) msg 138cc AJL/End 139c 140#include "bas_exndcf_sfn.fh" 141#include "bas_ibs_sfn.fh" 142c 143c check initialization and shells 144c 145 if (.not.int_chk_init('int_1eov')) 146 & call errquit('int_1eov: int_init was not called' ,0, 147 & INT_ERR) 148c 149 shells_ok = int_chk_sh(i_basis,ish) 150 shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh) 151 if (.not.shells_ok) 152 & call errquit('int_1eov: invalid contraction/shell',0, 153 & BASIS_ERR) 154c 155 ibas = i_basis + BASIS_HANDLE_OFFSET 156 jbas = j_basis + BASIS_HANDLE_OFFSET 157 oNR = .true. 158 oirel = .false. 159 ojrel = .false. 160 orel = .false. 161 canAB = .false. 162c 163 if (dyall_mod_dir) then 164c 165c get basis set handles; relativistic integral option only valid 166c if both ibas and jbas are the ao basis. 167c 168 lbas = lc_bsh + BASIS_HANDLE_OFFSET 169 sbas = sc_bsh + BASIS_HANDLE_OFFSET 170 abas = ao_bsh + BASIS_HANDLE_OFFSET 171 orel = ibas .eq. abas .and. jbas .eq. abas 172 end if 173c 174c i shell 175c 176 ucont = (sf_ibs_cn2ucn(ish,ibas)) 177c 178c check for relativistic shell 179c 180 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then 181 oirel = .true. 182 isbas = sbas 183 uconts = ao_to_ls(ucont) 184 if (uconts .eq. 0) call errquit ( 185 & 'int_1estv: no relativistic pointer',911, INT_ERR) 186 if (nesc_1e_approx) then 187 ibas = lbas 188 ucont = uconts 189 end if 190 else 191 uconts = ucont 192 isbas = ibas 193 end if 194c 195 Li = infbs_cont(CONT_TYPE ,ucont,ibas) 196 i_prim = infbs_cont(CONT_NPRIM,ucont,ibas) 197 i_gen = infbs_cont(CONT_NGEN ,ucont,ibas) 198 i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 199 i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas) 200 i_cent = (sf_ibs_cn2ce(ish,ibas)) 201 i_geom = ibs_geom(ibas) 202 i_icfpS = infbs_cont(CONT_ICFP ,uconts,isbas) 203c 204c j shell 205c 206 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 207c 208c check for relativistic shell 209c 210 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then 211 ojrel = .true. 212 jsbas = sbas 213 uconts = ao_to_ls(ucont) 214 if (uconts .eq. 0) call errquit ( 215 & 'int_1estv: no relativistic pointer',911, INT_ERR) 216 if (nesc_1e_approx) then 217 jbas = lbas 218 ucont = uconts 219 end if 220 else 221 uconts = ucont 222 jsbas = jbas 223 end if 224c 225 Lj = infbs_cont(CONT_TYPE ,ucont,jbas) 226 j_prim = infbs_cont(CONT_NPRIM,ucont,jbas) 227 j_gen = infbs_cont(CONT_NGEN ,ucont,jbas) 228 j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas) 229 j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas) 230 j_cent = (sf_ibs_cn2ce(jsh,jbas)) 231 j_geom = ibs_geom(jbas) 232 j_icfpS = infbs_cont(CONT_ICFP ,uconts,jsbas) 233c 234 oNR = .not.(oirel.and.ojrel) 235 orel = oirel.or.ojrel 236c 237 if (i_geom.ne.j_geom.and.WarnP.eq.0) then 238 write(luout,*) 239 & 'int_1eov: WARNING: possible geometry inconsistency' 240 write(luout,*)'i_basis geometry handle:',i_geom 241 write(luout,*)'j_basis geometry handle:',j_geom 242 WarnP = 1 243 endif 244 245 ohnd_ok = cando_hnd_1e(i_basis,ish,0) 246 & .and. cando_hnd_1e(j_basis,jsh,0) 247 & .and. (.not.geom_any_finuc (i_geom)) 248 & .and. (.not.geom_any_finuc (j_geom)) 249 onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0) 250 onw1e_ok = cando_nw_1e(i_basis,ish,0) 251 & .and. cando_nw_1e(j_basis,jsh,0) 252 if (orel) then 253 call rel_onel ( 254 & coords(1,i_cent,i_geom), 255 & dbl_mb(mb_exndcf(i_iexp,ibas)), 256 & dbl_mb(mb_exndcf(i_icfp,ibas)), 257 & dbl_mb(mb_exndcf(i_icfpS,isbas)),i_prim,i_gen,Li, 258 & coords(1,j_cent,j_geom), 259 & dbl_mb(mb_exndcf(j_iexp,jbas)), 260 & dbl_mb(mb_exndcf(j_icfp,jbas)), 261 & dbl_mb(mb_exndcf(j_icfpS,jsbas)),j_prim,j_gen,Lj, 262 & coords(1,1,i_geom),charge(1,i_geom), 263 & geom_invnucexp(1,i_geom),ncenter(i_geom), 264 & S,T,V,lstv,doS,doT,doV,canAB,onw_ok,ohnd_ok,oNR,.false., 265 & scr,lscr,rel_dbg,rel_typ) 266 else if (ohnd_ok) then 267 call hnd_stvint( 268 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 269 & dbl_mb(mb_exndcf(i_icfp,ibas)), 270 & i_prim, i_gen, Li, 271 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 272 & dbl_mb(mb_exndcf(j_icfp,jbas)), 273 & j_prim, j_gen, Lj, 274 & coords(1,1,i_geom),charge(1,i_geom),ncenter(i_geom), 275 & S,T,V,lstv,doS,doT,doV,scr,lscr) 276c 277 elseif (onw1e_ok) then 278 call int_hf1sp( 279 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 280 & dbl_mb(mb_exndcf(i_icfp,ibas)), 281 & i_prim, i_gen, Li, i_cent, 282 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 283 & dbl_mb(mb_exndcf(j_icfp,jbas)), 284 & j_prim, j_gen, Lj, j_cent, 285 & coords(1,1,i_geom),charge(1,i_geom), 286 & geom_invnucexp(1,i_geom),ncenter(i_geom), 287 & S,T,V,lstv,doS,doT,doV,canAB,.false., 288cc AJL/Begin 289c & scr,lscr,'int_1eov') 290 & scr,lscr,msg) 291cc AJL/End 292 elseif (onw_ok) then 293 call hf1( 294 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 295 & dbl_mb(mb_exndcf(i_icfp,ibas)), i_prim, i_gen, Li, 296 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 297 & dbl_mb(mb_exndcf(j_icfp,jbas)), j_prim, j_gen, Lj, 298 & coords(1,1,i_geom),charge(1,i_geom), 299 & geom_invnucexp(1,i_geom),ncenter(i_geom), 300 & S,T,V,lstv,doS,doT,doV,canAB,.false., 301 & scr,lscr) 302 else 303 call errquit('int_1eov: could not do hnd, sp or nw integrals', 304 & 0, INT_ERR) 305 endif 306c 307* We now have the cartesian integral block(s) (jlo:jhi,ilo:ihi) 308* 309 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 310 if (.not.any_spherical) return 311c 312c ... reset general contractions for sp shells to 1 since they are handled 313c as a block of 4. Since int_nbf_* arrays are set to the appropriate size. 314c 315 if (li.eq.-1) i_gen = 1 316 if (lj.eq.-1) j_gen = 1 317c 318 if (bas_spherical(ibas).and.bas_spherical(jbas)) then 319*... transform both i and j integrals 320 i_nbf_x = int_nbf_x(Li) 321 i_nbf_s = int_nbf_s(Li) 322 j_nbf_x = int_nbf_x(Lj) 323 j_nbf_s = int_nbf_s(Lj) 324c 325 if (doS) call spcart_tran1e(S,scr, 326 & j_nbf_x,i_nbf_x,Lj,j_gen, 327 & j_nbf_s,i_nbf_s,Li,i_gen, 328 & .false.) 329 if (doT) call spcart_tran1e(T,scr, 330 & j_nbf_x,i_nbf_x,Lj,j_gen, 331 & j_nbf_s,i_nbf_s,Li,i_gen, 332 & .false.) 333 if (doV) call spcart_tran1e(V,scr, 334 & j_nbf_x,i_nbf_x,Lj,j_gen, 335 & j_nbf_s,i_nbf_s,Li,i_gen, 336 & .false.) 337 else if (bas_spherical(ibas)) then 338*.. transform on i component 339 i_nbf_x = int_nbf_x(Li) 340 i_nbf_s = int_nbf_s(Li) 341 j_nbf_x = int_nbf_x(Lj) 342 j_nbf_s = j_nbf_x 343 if (doS) call spcart_tran1e(S,scr, 344 & j_nbf_x,i_nbf_x,0,j_gen, 345 & j_nbf_s,i_nbf_s,Li,i_gen, 346 & .false.) 347 if (doT) call spcart_tran1e(T,scr, 348 & j_nbf_x,i_nbf_x,0,j_gen, 349 & j_nbf_s,i_nbf_s,Li,i_gen, 350 & .false.) 351 if (doV) call spcart_tran1e(V,scr, 352 & j_nbf_x,i_nbf_x,0,j_gen, 353 & j_nbf_s,i_nbf_s,Li,i_gen, 354 & .false.) 355 else if (bas_spherical(jbas)) then 356*.. transform on j component 357 i_nbf_x = int_nbf_x(Li) 358 i_nbf_s = i_nbf_x 359 j_nbf_x = int_nbf_x(Lj) 360 j_nbf_s = int_nbf_s(Lj) 361 if (doS) call spcart_tran1e(S,scr, 362 & j_nbf_x,i_nbf_x,Lj,j_gen, 363 & j_nbf_s,i_nbf_s,0,i_gen, 364 & .false.) 365 if (doT) call spcart_tran1e(T,scr, 366 & j_nbf_x,i_nbf_x,Lj,j_gen, 367 & j_nbf_s,i_nbf_s,0,i_gen, 368 & .false.) 369 if (doV) call spcart_tran1e(V,scr, 370 & j_nbf_x,i_nbf_x,Lj,j_gen, 371 & j_nbf_s,i_nbf_s,0,i_gen, 372 & .false.) 373 else 374 call errquit( 375 & 'int_1eov: should never reach transform blocked else',911, 376 & INT_ERR) 377 endif 378 return 379c 380 end 381C> @} 382c---------------------------------------------------------------------- 383c 384c Wrapper routines to compute different classes of integrals. 385c All of these call int_1estv. 386c 387c---------------------------------------------------------------------- 388* 389C> \ingroup nwint 390C> @{ 391C> \brief Compute the overlap integrals 392C> 393C> Compute the overlap integrals between two shells. These integrals 394C> are defined as: 395C> \f{eqnarray*}{ 396C> S & = & ({\mu}|{\nu}) \\\\ 397C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} 398C> \f} 399c:tex-% this is part of the API Standard Integral routines. 400c:tex-\subsection{int\_1eov} 401c:tex-This routine computes the 1 electron overlap integrals ($S$): 402c:tex-\begin{eqnarray*} 403c:tex-S & = & ({\mu}|{\nu}) \\ 404c:tex- & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} 405c:tex-\end{eqnarray*} 406c:tex- 407c:tex-{\it Syntax:} 408c:tex-\begin{verbatim} 409 subroutine int_1eov(i_basis,ish,j_basis,jsh,lscr,scr,lov,Ov) 410#include "util.fh" 411c:tex-\end{verbatim} 412c 413c::passed 414c:tex-\begin{verbatim} 415 integer i_basis !< [Input] basis set handle for ish 416 integer ish !< [Input] i shell/contraction 417 integer j_basis !< [Input] basis set handle for jsh 418 integer jsh !< [Input] j shell/contraction 419 integer lscr !< [Input] length of scratch array 420 double precision scr(lscr) !< [Scratch] scratch array 421 integer lov !< [Input] length of Ov buffer 422 double precision Ov(lov) !< [Output] overlap integrals 423c:tex-\end{verbatim} 424c 425 call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lov,Ov,scr,scr, 426cc AJL/Begin 427c & .true.,.false.,.false.) 428 & .true.,.false.,.false.,'int_1eov') 429c doS doT doV msg 430cc AJL/End 431 return 432 end 433c---------------------------------------------------------------------- 434* 435C> \brief Calculates the 1-electron kinetic energy integrals 436C> 437C> This routine calculates the 1-electron kinetic energy integrals 438C> between a pair of specified shells. The integrals are defined as: 439C> \f{eqnarray*}{ 440C> T & = & ({\mu}|\frac{-1}{2}\nabla^{2}|{\nu}) \\\\ 441C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\frac{-1}{2}\nabla^{2}(r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} 442C> \f} 443c:tex-% this is part of the API Standard Integral routines. 444c:tex-\subsection{int\_1eke} 445c:tex-This routine computes the 1 electron kinetic energy integrals, ($T$).: 446c:tex-\begin{eqnarray*} 447c:tex-T & = & ({\mu}|\frac{-1}{2}\nabla^{2}|{\nu}) \\ 448c:tex- & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\frac{-1}{2}\nabla^{2}(r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} 449c:tex-\end{eqnarray*} 450c:tex- 451c:tex-{\it Syntax:} 452c:tex-\begin{verbatim} 453 subroutine int_1eke(i_basis,ish,j_basis,jsh,lscr,scr,lke,Ke) 454#include "util.fh" 455c:tex-\end{verbatim} 456c 457c::passed 458c:tex-\begin{verbatim} 459 integer i_basis !< [Input] basis set handle for ish 460 integer ish !< [Input] i shell/contraction 461 integer j_basis !< [Input] basis set handle for jsh 462 integer jsh !< [Input] j shell/contraction 463 integer lscr !< [Input] length of scratch array 464 double precision scr(lscr) !< [Scratch] scratch array 465 integer lke !< [Input] length of Ke buffer 466 double precision Ke(lke) !< [Output] kinetic energy integrals 467c:tex-\end{verbatim} 468c 469 call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lke,scr,Ke,scr, 470cc AJL/Begin 471c & .false.,.true.,.false.) 472 & .false.,.true.,.false.,'int_1eke') 473c doS doT doV msg 474cc AJL/End 475 return 476 end 477c---------------------------------------------------------------------- 478* 479C> \brief Computes the nuclear attraction integrals 480C> 481C> This routine calculates the 1-electron nuclear attraction integrals 482C> for a pair of shells. If any ECPs are present then the ECP terms 483C> are directly summed into these integrals. 484C> 485C> The integrals are defined as: 486C> \f{eqnarray*}{ 487C> V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\\\ 488C> & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 489C> {-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 490C> \f} 491c:tex-% this is part of the API Standard Integral routines. 492c:tex-\subsection{int\_1epe} 493c:tex-This routine computes the 1 electron potential integrals, ($V$): 494c:tex-\begin{eqnarray*} 495c:tex-V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\ 496c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 497c:tex-{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 498c:tex-\end{eqnarray*} 499c:tex-If an ECP is defined then the ECP integral contributions are summed 500c:tex-directly into the V integrals. 501c:tex- 502c:tex-{\it Syntax:} 503c:tex-\begin{verbatim} 504 subroutine int_1epe(i_basis,ish,j_basis,jsh,lscr,scr,lpe,Pe) 505#include "util.fh" 506c:tex-\end{verbatim} 507c 508c::passed 509c:tex-\begin{verbatim} 510 integer i_basis !< [Input] basis set handle for ish 511 integer ish !< [Input] i shell/contraction 512 integer j_basis !< [Input] basis set handle for jsh 513 integer jsh !< [Input] j shell/contraction 514 integer lscr !< [Input] length of scratch array 515 double precision scr(lscr) !< [Scratch] scratch array 516 integer lpe !< [Input] length of Pe buffer 517 double precision Pe(lpe) !< [Output] kinetic energy integrals 518c:tex-\end{verbatim} 519c 520 call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr,scr,Pe, 521cc AJL/Begin 522c & .false.,.false.,.true.) 523 & .false.,.false.,.true.,'int_1epe') 524c doS doT doV msg 525cc AJL/End 526 return 527 end 528C> @} 529cc AJL/Begin/SPIN ECPs 530c---------------------------------------------------------------------- 531* 532C> \brief Computes the nuclear attraction integrals 533C> 534C> This routine calculates the 1-electron nuclear attraction integrals 535C> for a pair of shells. If any ECPs are present then the ECP terms 536C> are directly summed into these integrals. 537C> 538C> The integrals are defined as: 539C> \f{eqnarray*}{ 540C> V & = & 541C({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\\\ 542C> & = & \int_{-\infty}^{\infty} 543Cg_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 544C> {-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 545C> \f} 546c:tex-% this is part of the API Standard Integral routines. 547c:tex-\subsection{int\_1epe} 548c:tex-This routine computes the 1 electron potential integrals, ($V$): 549c:tex-\begin{eqnarray*} 550c:tex-V & = &({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\ 551c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 552c:tex-{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 553c:tex-\end{eqnarray*} 554c:tex-If a *BETA* ECP is defined then the ECP integral contributions are summed 555c:tex-directly into the V integrals. 556c:tex- 557c:tex-{\it Syntax:} 558c:tex-\begin{verbatim} 559 subroutine int_1epe_beta(i_basis,ish,j_basis,jsh,lscr,scr,lpe,Pe) 560#include "util.fh" 561c:tex-\end{verbatim} 562c 563c::passed 564c:tex-\begin{verbatim} 565 integer i_basis !< [Input] basis set handle for ish 566 integer ish !< [Input] i shell/contraction 567 integer j_basis !< [Input] basis set handle for jsh 568 integer jsh !< [Input] j shell/contraction 569 integer lscr !< [Input] length of scratch array 570 double precision scr(lscr) !< [Scratch] scratch array 571 integer lpe !< [Input] length of Pe buffer 572 double precision Pe(lpe) !< [Output] kinetic energy integrals 573c:tex-\end{verbatim} 574c 575 call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr,scr,Pe, 576 & .false.,.false.,.true.,'int_1epe_beta') 577c doS doT doV msg 578 return 579 end 580cc AJL/End 581C> @} 582c---------------------------------------------------------------------- 583* 584C> \ingroup nwint 585C> @{ 586C> \brief Compute the 1-electron Hamiltonian 587C> 588C> This routine computes the 1 electron hamiltonian \f$H1\f$. 589C> \f{eqnarray*}{ 590C> H1 & = & T + V \\\\ 591C> T & = & ({\mu}|\frac{-1}{2}\nabla^{2}|{\nu}) \\\\ 592C> & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\frac{-1}{2} 593C> \nabla^{2}(r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\\\ 594C> V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\\\ 595C> & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 596C> {-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 597C> \f} 598C> 599C> If an ECP is defined then the ECP integral contributions are summed 600C> directly into the \f$H1\f$ integrals. 601C> 602C> If a relativistic basis is defined then the one-electron integrals for 603C> the case where both shells are relativistic are the modified integrals. 604C> 605c:tex-% this is part of the API Standard Integral routines. 606c:tex-\subsection{int\_1eh1} 607c:tex-This routine computes the 1 electron hamiltonian, ($H1$). 608c:tex-\begin{eqnarray*} 609c:tex-H1 & = & T + V \\ 610c:tex-T & = & ({\mu}|\frac{-1}{2}\nabla^{2}|{\nu}) \\ 611c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\frac{-1}{2} 612c:tex-\nabla^{2}(r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\ 613c:tex-V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\ 614c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac 615c:tex-{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 616c:tex-\end{eqnarray*} 617c:tex- 618c:tex-If an ECP is defined then the ECP integral contributions are summed 619c:tex-directly into the $H1$ integrals. 620c:tex- 621c:tex-If a relativistic basis is defined then the one-electron integrals for 622c:tex-the case where both shells are relativistic are the modified integrals. 623c:tex- 624c:tex-{\it Syntax:} 625c:tex-\begin{verbatim} 626 subroutine int_1eh1(i_basis,ish,j_basis,jsh,lscr,scr,lh1,H1) 627#include "util.fh" 628c:tex-\end{verbatim} 629#include "apiP.fh" 630c 631c::functions 632 logical cando_hnd_1e 633 external cando_hnd_1e 634c::passed 635c:tex-\begin{verbatim} 636 integer i_basis !< [Input] basis set handle for ish 637 integer ish !< [Input] i shell/contraction 638 integer j_basis !< [Input] basis set handle for jsh 639 integer jsh !< [Input] j shell/contraction 640 integer lscr !< [Input] length of scratch array 641 double precision scr(lscr) !< [Scratch] scratch array 642 integer lh1 !< [Input] length of H1 buffer. 643 double precision H1(lh1) !< [Output] one electron 644c:tex-\end{verbatim} 645c 646 call int_1estv(i_basis,ish,j_basis,jsh,lscr-lh1,scr(lh1+1),lh1, 647cc AJL/Begin 648c & scr,scr,H1,.false.,.true.,.true.) 649 & scr,scr,H1,.false.,.true.,.true.,'int_1eh1') 650c doS doT doV msg 651cc AJL/End 652 do i = 1,lh1 653 H1(i) = H1(i)+scr(i) 654 end do 655 return 656 end 657c---------------------------------------------------------------------- 658* 659C> \brief Compute all 1-electron integrals 660C> 661C> This routine computes the 1 electron integrals \f$S\f$, \f$T\f$, and \f$V\f$: 662C> \f{eqnarray*}{ 663C> S & = & ({\mu}|{\nu}) \\\\ 664C> & = & \int_{{-}\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\\\ 665C> T & = & ({\mu}|-\frac{1}{2}{\nabla^{2}}|{\nu}) \\\\ 666C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\frac{-1}{2}{\nabla^{2}}(r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\\\ 667C> V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\\\ 668C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-Z_{\alpha}} 669C> {|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 670C> \f} 671C> 672c:tex- 673c:tex-% this is part of the API Standard Integral routines. 674c:tex-\subsection{int\_1eall} 675c:tex-This routine computes the 1 electron integrals S, T, and V: 676c:tex-\begin{eqnarray*} 677c:tex-S & = & ({\mu}|{\nu}) \\ 678c:tex- & = & \int_{{-}\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\ 679c:tex-T & = & ({\mu}|-\frac{1}{2}{\nabla^{2}}|{\nu}) \\ 680c:tex- & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\frac{-1}{2}{\nabla^{2}}(r_{1})g_{\nu}(X_{\nu},r_{1})dr_{1} \\ 681c:tex-V & = & ({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\ 682c:tex- & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-Z_{\alpha}} 683c:tex-{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 684c:tex-\end{eqnarray*} 685c:tex- 686c:tex- 687c:tex-{\it Syntax:} 688c:tex-\begin{verbatim} 689 subroutine int_1eall(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V) 690#include "util.fh" 691c:tex-\end{verbatim} 692c 693c::passed 694c:tex-\begin{verbatim} 695 integer i_basis !< [Input] basis set handle for ish 696 integer ish !< [Input] i shell/contraction 697 integer j_basis !< [Input] basis set handle for jsh 698 integer jsh !< [Input] j shell/contraction 699 integer lscr !< [Input] length of scratch array 700 double precision scr(lscr) !< [Scratch] scratch array 701 integer lstv !< [Input] length of one electron buffers 702 double precision T(lstv) !< [Output] kinetic integral buffer 703 double precision V(lstv) !< [Output] potential integral buffer 704 double precision S(lstv) !< [Output] overlap integral buffer 705c:tex-\end{verbatim} 706c 707 call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V, 708cc AJL/Begin 709c & .true.,.true.,.true.) 710 & .true.,.true.,.true.,'int_1eall') 711c doS doT doV msg 712cc AJL/End 713 return 714 end 715c---------------------------------------------------------------------- 716C> \brief Compute the 1-electron COSMO embedding potential 717C> 718C> This routine computes the 1 electron integrals \f$V\f$: 719C> \f{eqnarray*}{ 720C> V & = & ({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\ 721C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}} 722C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 723C> \f} 724C> where \f$ q_\alpha \f$ is a point charge situated at position \f$ r_\alpha \f$. 725C> The COSMO embedding point charges are accessed in int_1eefc. 726C> 727 subroutine int_1epot(i_basis,ish,j_basis,jsh,lscr,scr,lve,Ve) 728#include "util.fh" 729c 730 integer i_basis !< [Input] basis set handle for ish 731 integer ish !< [Input] i shell/contraction 732 integer j_basis !< [Input] basis set handle for jsh 733 integer jsh !< [Input] j shell/contraction 734 integer lscr !< [Input] length of scratch array 735 double precision scr(lscr) !< [Scratch] scratch array 736 integer lve !< [Input] length of Ve buffer 737 double precision Ve(lve) !< [Output] potential energy integrals 738c 739 call int_1eefc(i_basis,ish,j_basis,jsh,lscr,scr,lve,scr,scr,Ve, 740 & .false.,.false.,.true.) 741c doS doT doV 742 return 743 end 744c---------------------------------------------------------------------- 745C> \brief Compute the 1-electron integrals but the potential is due to 746C> COSMO embedding charges 747C> 748C> This routine computes the 1 electron integrals \f$V\f$: 749C> \f{eqnarray*}{ 750C> V & = & ({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\ 751C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}} 752C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 753C> \f} 754C> where \f$ q_\alpha \f$ is a point charge situated at position \f$ r_\alpha \f$. 755C> The location of the COSMO embedding point charges is specified in the common 756C> block include from prop.fh. 757C> 758 subroutine int_1eefc(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V, 759 & doS,doT,doV) 760 implicit none 761#include "nwc_const.fh" 762#include "errquit.fh" 763#include "basP.fh" 764#include "basdeclsP.fh" 765#include "geomP.fh" 766#include "geom.fh" 767#include "bq.fh" 768#include "geobasmapP.fh" 769#include "mafdecls.fh" 770#include "bas_exndcf_dec.fh" 771#include "bas_ibs_dec.fh" 772#include "int_nbf.fh" 773#include "stdio.fh" 774#include "apiP.fh" 775#include "rel_nwc.fh" 776#include "util.fh" 777C 778C The following include block contains the geometry handle for the 779C geometry and charges of the solvent for COSMO. This is a clear violation 780C of the integral API, but for now the cleanest solution 781C 782#include "prop.fh" 783c::external subroutines used 784c... errquit 785c::functions 786 logical cando_hnd_1e 787 logical cando_nw_1e 788 logical cando_nw 789 logical int_chk_init 790 logical int_chk_sh 791 external int_chk_init 792 external int_chk_sh 793 external cando_hnd_1e 794 external cando_nw_1e 795 external cando_nw 796c::passed 797c:tex-\begin{verbatim} 798 integer i_basis !< [Input] basis set handle for ish 799 integer ish !< [Input] i shell/contraction 800 integer j_basis !< [Input] basis set handle for jsh 801 integer jsh !< [Input] j shell/contraction 802 integer lscr !< [Input] length of scratch array 803 integer lstv !< [Input] length of integral buffer 804 double precision scr(lscr) !< [Scratch] scratch array 805 double precision S(lstv) !< [Output] overlap integrals 806 double precision T(lstv) !< [Output] kinetic energy integrals 807 double precision V(lstv) !< [Output] potential energy integrals 808 logical doS !< [Input] flag for overlap integrals 809 logical doT !< [Input] flag for kinetic energy integrals 810 logical doV !< [Input] flag for potential energy integrals 811c:tex-\end{verbatim} 812c::local 813 logical ohnd_ok, onw_ok, onw1e_ok 814 logical shells_ok, orel, oirel, ojrel, oNR, canAB 815 integer i_geom, j_geom, ibas, jbas, ucont, uconts 816 integer lbas, sbas, abas, isbas, jsbas 817 integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_icfpS 818 integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_icfpS 819 integer iefc_c, iefc_q, nefc 820 integer k_efc_inv 821c 822 logical any_spherical 823 integer i_nbf_x, j_nbf_x 824 integer i_nbf_s, j_nbf_s 825c 826 integer WarnP 827 save WarnP 828 data WarnP /0/ 829c 830#include "bas_exndcf_sfn.fh" 831#include "bas_ibs_sfn.fh" 832c 833c check initialization and shells 834c 835 if (.not.int_chk_init('int_1eefc')) 836 & call errquit('int_1eefc: int_init was not called' ,0, 837 & INT_ERR) 838c 839 shells_ok = int_chk_sh(i_basis,ish) 840 shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh) 841 if (.not.shells_ok) 842 & call errquit('int_1eefc: invalid contraction/shell',0, 843 & BASIS_ERR) 844c 845 ibas = i_basis + BASIS_HANDLE_OFFSET 846 jbas = j_basis + BASIS_HANDLE_OFFSET 847 oNR = .true. 848 oirel = .false. 849 ojrel = .false. 850 orel = .false. 851 canAB = .false. 852c 853 if (dyall_mod_dir) then 854c 855c get basis set handles; relativistic integral option only valid 856c if both ibas and jbas are the ao basis. 857c 858 lbas = lc_bsh + BASIS_HANDLE_OFFSET 859 sbas = sc_bsh + BASIS_HANDLE_OFFSET 860 abas = ao_bsh + BASIS_HANDLE_OFFSET 861 orel = ibas .eq. abas .and. jbas .eq. abas 862 end if 863c 864c i shell 865c 866 ucont = (sf_ibs_cn2ucn(ish,ibas)) 867c 868c check for relativistic shell 869c 870 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then 871 oirel = .true. 872 isbas = sbas 873 uconts = ao_to_ls(ucont) 874 if (uconts .eq. 0) call errquit ( 875 & 'int_1eefc: no relativistic pointer',911, INT_ERR) 876 if (nesc_1e_approx) then 877 ibas = lbas 878 ucont = uconts 879 end if 880 else 881 uconts = ucont 882 isbas = ibas 883 end if 884c 885 Li = infbs_cont(CONT_TYPE ,ucont,ibas) 886 i_prim = infbs_cont(CONT_NPRIM,ucont,ibas) 887 i_gen = infbs_cont(CONT_NGEN ,ucont,ibas) 888 i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 889 i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas) 890 i_cent = (sf_ibs_cn2ce(ish,ibas)) 891 i_geom = ibs_geom(ibas) 892 i_icfpS = infbs_cont(CONT_ICFP ,uconts,isbas) 893c 894c j shell 895c 896 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 897c 898c check for relativistic shell 899c 900 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then 901 ojrel = .true. 902 jsbas = sbas 903 uconts = ao_to_ls(ucont) 904 if (uconts .eq. 0) call errquit ( 905 & 'int_1eefc: no relativistic pointer',911, INT_ERR) 906 if (nesc_1e_approx) then 907 jbas = lbas 908 ucont = uconts 909 end if 910 else 911 uconts = ucont 912 jsbas = jbas 913 end if 914c 915 Lj = infbs_cont(CONT_TYPE ,ucont,jbas) 916 j_prim = infbs_cont(CONT_NPRIM,ucont,jbas) 917 j_gen = infbs_cont(CONT_NGEN ,ucont,jbas) 918 j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas) 919 j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas) 920 j_cent = (sf_ibs_cn2ce(jsh,jbas)) 921 j_geom = ibs_geom(jbas) 922 j_icfpS = infbs_cont(CONT_ICFP ,uconts,jsbas) 923c 924 oNR = .not.(oirel.and.ojrel) 925 orel = oirel.or.ojrel 926c 927 if (i_geom.ne.j_geom.and.WarnP.eq.0) then 928 write(luout,*) 929 & 'int_1eefc: WARNING: possible geometry inconsistency' 930 write(luout,*)'i_basis geometry handle:',i_geom 931 write(luout,*)'j_basis geometry handle:',j_geom 932 WarnP = 1 933 endif 934c 935c ----- current cosmo restriction ... ----- 936c 937 if(i_geom.ne.j_geom) then 938 write(luout,*) 'i_geom andr j_geom is noteq' 939 call errquit('int_1eefc: i_geom and j_geom must be the same',0, 940 & INT_ERR) 941 endif 942c 943 if (.not.bq_ncenter(cosmo_bq_efc,nefc)) 944 + call errquit("int_1eefc: bq_ncenter failed",nefc,uerr) 945 if (.not.bq_index_coord(cosmo_bq_efc,iefc_c)) 946 + call errquit("int_1eefc: bq_index_coord failed",iefc_c,uerr) 947 if (.not.bq_index_charge(cosmo_bq_efc,iefc_q)) 948 + call errquit("int_1eefc: bq_index_charge failed",iefc_q,uerr) 949c 950c Pick up the array we created and zeroed in cosmo_initialize 951c 952 if (.not.bq_index_charge(cosmo_bq_invnuc,k_efc_inv)) 953 + call errquit("int_1eefc: bq_index_charge failed",iefc_q,uerr) 954c 955 ohnd_ok = cando_hnd_1e(i_basis,ish,0) 956 & .and. cando_hnd_1e(j_basis,jsh,0) 957 & .and. (.not.geom_any_finuc (i_geom)) 958 & .and. (.not.geom_any_finuc (j_geom)) 959 onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0) 960 onw1e_ok = cando_nw_1e(i_basis,ish,0) 961 & .and. cando_nw_1e(j_basis,jsh,0) 962c 963 if (orel) then 964 call rel_onel ( 965 & coords(1,i_cent,i_geom), 966 & dbl_mb(mb_exndcf(i_iexp,ibas)), 967 & dbl_mb(mb_exndcf(i_icfp,ibas)), 968 & dbl_mb(mb_exndcf(i_icfpS,isbas)),i_prim,i_gen,Li, 969 & coords(1,j_cent,j_geom), 970 & dbl_mb(mb_exndcf(j_iexp,jbas)), 971 & dbl_mb(mb_exndcf(j_icfp,jbas)), 972 & dbl_mb(mb_exndcf(j_icfpS,jsbas)),j_prim,j_gen,Lj, 973 & dbl_mb(iefc_c),dbl_mb(iefc_q), 974 & dbl_mb(k_efc_inv),nefc, 975 & S,T,V,lstv,doS,doT,doV,canAB,onw_ok,ohnd_ok,oNR,.false., 976 & scr,lscr,rel_dbg,rel_typ) 977 else if (ohnd_ok) then 978 call hnd_stvint( 979 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 980 & dbl_mb(mb_exndcf(i_icfp,ibas)), 981 & i_prim, i_gen, Li, 982 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 983 & dbl_mb(mb_exndcf(j_icfp,jbas)), 984 & j_prim, j_gen, Lj, 985 & dbl_mb(iefc_c),dbl_mb(iefc_q),nefc, 986 & S,T,V,lstv,doS,doT,doV,scr,lscr) 987c 988 elseif (onw1e_ok) then 989 call int_hf1sp( 990 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 991 & dbl_mb(mb_exndcf(i_icfp,ibas)), 992 & i_prim, i_gen, Li, i_cent, 993 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 994 & dbl_mb(mb_exndcf(j_icfp,jbas)), 995 & j_prim, j_gen, Lj, j_cent, 996 & dbl_mb(iefc_c),dbl_mb(iefc_q), 997 & dbl_mb(k_efc_inv),nefc, 998 & S,T,V,lstv,doS,doT,doV,canAB,.false., 999 & scr,lscr,'int_1eefc') 1000 elseif (onw_ok) then 1001 call hf1( 1002 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1003 & dbl_mb(mb_exndcf(i_icfp,ibas)), i_prim, i_gen, Li, 1004 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1005 & dbl_mb(mb_exndcf(j_icfp,jbas)), j_prim, j_gen, Lj, 1006 & dbl_mb(iefc_c),dbl_mb(iefc_q), 1007 & dbl_mb(k_efc_inv),nefc, 1008 & S,T,V,lstv,doS,doT,doV,canAB,.false., 1009 & scr,lscr) 1010 else 1011 call errquit('int_1eefc: could not do hnd, sp or nw integrals', 1012 & 0, INT_ERR) 1013 endif 1014c 1015* We now have the cartesian integral block(s) (jlo:jhi,ilo:ihi) 1016* 1017 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 1018 if (.not.any_spherical) return 1019c 1020c ... reset general contractions for sp shells to 1 since they are handled 1021c as a block of 4. Since int_nbf_* arrays are set to the appropriate size. 1022c 1023 if (li.eq.-1) i_gen = 1 1024 if (lj.eq.-1) j_gen = 1 1025c 1026 if (bas_spherical(ibas).and.bas_spherical(jbas)) then 1027*... transform both i and j integrals 1028 i_nbf_x = int_nbf_x(Li) 1029 i_nbf_s = int_nbf_s(Li) 1030 j_nbf_x = int_nbf_x(Lj) 1031 j_nbf_s = int_nbf_s(Lj) 1032c 1033 if (doS) call spcart_tran1e(S,scr, 1034 & j_nbf_x,i_nbf_x,Lj,j_gen, 1035 & j_nbf_s,i_nbf_s,Li,i_gen, 1036 & .false.) 1037 if (doT) call spcart_tran1e(T,scr, 1038 & j_nbf_x,i_nbf_x,Lj,j_gen, 1039 & j_nbf_s,i_nbf_s,Li,i_gen, 1040 & .false.) 1041 if (doV) call spcart_tran1e(V,scr, 1042 & j_nbf_x,i_nbf_x,Lj,j_gen, 1043 & j_nbf_s,i_nbf_s,Li,i_gen, 1044 & .false.) 1045 else if (bas_spherical(ibas)) then 1046*.. transform on i component 1047 i_nbf_x = int_nbf_x(Li) 1048 i_nbf_s = int_nbf_s(Li) 1049 j_nbf_x = int_nbf_x(Lj) 1050 j_nbf_s = j_nbf_x 1051 if (doS) call spcart_tran1e(S,scr, 1052 & j_nbf_x,i_nbf_x,0,j_gen, 1053 & j_nbf_s,i_nbf_s,Li,i_gen, 1054 & .false.) 1055 if (doT) call spcart_tran1e(T,scr, 1056 & j_nbf_x,i_nbf_x,0,j_gen, 1057 & j_nbf_s,i_nbf_s,Li,i_gen, 1058 & .false.) 1059 if (doV) call spcart_tran1e(V,scr, 1060 & j_nbf_x,i_nbf_x,0,j_gen, 1061 & j_nbf_s,i_nbf_s,Li,i_gen, 1062 & .false.) 1063 else if (bas_spherical(jbas)) then 1064*.. transform on j component 1065 i_nbf_x = int_nbf_x(Li) 1066 i_nbf_s = i_nbf_x 1067 j_nbf_x = int_nbf_x(Lj) 1068 j_nbf_s = int_nbf_s(Lj) 1069 if (doS) call spcart_tran1e(S,scr, 1070 & j_nbf_x,i_nbf_x,Lj,j_gen, 1071 & j_nbf_s,i_nbf_s,0,i_gen, 1072 & .false.) 1073 if (doT) call spcart_tran1e(T,scr, 1074 & j_nbf_x,i_nbf_x,Lj,j_gen, 1075 & j_nbf_s,i_nbf_s,0,i_gen, 1076 & .false.) 1077 if (doV) call spcart_tran1e(V,scr, 1078 & j_nbf_x,i_nbf_x,Lj,j_gen, 1079 & j_nbf_s,i_nbf_s,0,i_gen, 1080 & .false.) 1081 else 1082 call errquit( 1083 & 'int_1eefc: should never reach transform blocked else',911, 1084 & INT_ERR) 1085 endif 1086 return 1087c 1088 end 1089 1090c---------------------------------------------------------------------- 1091C> \brief Compute the 1-electron Bq embedding potential 1092C> 1093C> This routine computes the 1 electron integrals \f$V\f$: 1094C> \f{eqnarray*}{ 1095C> V & = & ({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\ 1096C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}} 1097C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 1098C> \f} 1099C> where \f$ q_\alpha \f$ is a point charge situated at position \f$ r_\alpha \f$. 1100C> The Bq embedding point charges are accessed in int_1eefc1. 1101C> 1102 subroutine int_1epot1(i_basis,ish,j_basis,jsh,lscr,scr,lpe,Pe) 1103#include "util.fh" 1104c 1105 integer i_basis !< [Input] basis set handle for ish 1106 integer ish !< [Input] i shell/contraction 1107 integer j_basis !< [Input] basis set handle for jsh 1108 integer jsh !< [Input] j shell/contraction 1109 integer lscr !< [Input] length of scratch array 1110 double precision scr(lscr) !< [Scratch] scratch array 1111 integer lpe !< [Input] length of Pe buffer 1112 double precision Pe(lpe) !< [Output] potential energy integrals 1113c 1114 call int_1eefc1(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr,scr,Pe, 1115 & .false.,.false.,.true.) 1116c doS doT doV 1117 return 1118 end 1119c---------------------------------------------------------------------- 1120C> \brief Compute the 1-electron integrals but the potential is due to 1121C> Bq embedding charges 1122C> 1123C> This routine computes the 1 electron integrals \f$V\f$: 1124C> \f{eqnarray*}{ 1125C> V & = & ({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\ 1126C> & = & \int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}} 1127C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 1128C> \f} 1129C> where \f$ q_\alpha \f$ is a point charge situated at position \f$ r_\alpha \f$. 1130C> The location of the Bq embedding point charges is specified in a 1131C> Bq instance that has been activated elsewhere. 1132C> 1133 subroutine int_1eefc1(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V, 1134 & doS,doT,doV) 1135 implicit none 1136#include "nwc_const.fh" 1137#include "errquit.fh" 1138#include "basP.fh" 1139#include "basdeclsP.fh" 1140#include "geomP.fh" 1141#include "geom.fh" 1142#include "geobasmapP.fh" 1143#include "mafdecls.fh" 1144#include "bas_exndcf_dec.fh" 1145#include "bas_ibs_dec.fh" 1146#include "int_nbf.fh" 1147#include "stdio.fh" 1148#include "apiP.fh" 1149#include "rel_nwc.fh" 1150#include "util.fh" 1151C 1152C The following include block contains the geometry handle for the 1153C geometry and charges of the solvent for COSMO. This is a clear violation 1154C of the integral API, but for now the cleanest solution 1155C 1156#include "prop.fh" 1157c::external subroutines used 1158c... errquit 1159c::functions 1160 logical cando_hnd_1e 1161 logical cando_nw_1e 1162 logical cando_nw 1163 logical int_chk_init 1164 logical int_chk_sh 1165 external int_chk_init 1166 external int_chk_sh 1167 external cando_hnd_1e 1168 external cando_nw_1e 1169 external cando_nw 1170c::passed 1171c:tex-\begin{verbatim} 1172 integer i_basis !< [Input] basis set handle for ish 1173 integer ish !< [Input] i shell/contraction 1174 integer j_basis !< [Input] basis set handle for jsh 1175 integer jsh !< [Input] j shell/contraction 1176 integer lscr !< [Input] length of scratch array 1177 integer lstv !< [Input] length of integral buffer 1178 double precision scr(lscr) !< [Scratch] scratch array 1179 double precision S(lstv) !< [Output] overlap integrals 1180 double precision T(lstv) !< [Output] kinetic energy integrals 1181 double precision V(lstv) !< [Output] potential energy integrals 1182 logical doS !< [Input] flag for overlap integrals 1183 logical doT !< [Input] flag for kinetic energy integrals 1184 logical doV !< [Input] flag for potential energy integrals 1185c:tex-\end{verbatim} 1186c::local 1187 logical ohnd_ok, onw_ok, onw1e_ok 1188 logical shells_ok, orel, oirel, ojrel, oNR, canAB 1189 integer i_geom, j_geom, ibas, jbas, ucont, uconts 1190 integer lbas, sbas, abas, isbas, jsbas 1191 integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_icfpS 1192 integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_icfpS 1193c 1194 logical any_spherical 1195 integer i_nbf_x, j_nbf_x 1196 integer i_nbf_s, j_nbf_s 1197c 1198 integer WarnP 1199 save WarnP 1200 data WarnP /0/ 1201c 1202 integer bq_ncent 1203 integer i_qbq,i_cbq,h_xbq0,i_xbq0 1204c 1205#include "bas_exndcf_sfn.fh" 1206#include "bas_ibs_sfn.fh" 1207c 1208c check initialization and shells 1209c 1210 if (.not.int_chk_init('int_1eefc1')) 1211 & call errquit('int_1eefc1: int_init was not called' ,0, 1212 & INT_ERR) 1213c 1214 shells_ok = int_chk_sh(i_basis,ish) 1215 shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh) 1216 if (.not.shells_ok) 1217 & call errquit('int_1eefc1: invalid contraction/shell',0, 1218 & BASIS_ERR) 1219c 1220 ibas = i_basis + BASIS_HANDLE_OFFSET 1221 jbas = j_basis + BASIS_HANDLE_OFFSET 1222 oNR = .true. 1223 oirel = .false. 1224 ojrel = .false. 1225 orel = .false. 1226 canAB = .false. 1227c 1228 if (dyall_mod_dir) then 1229c 1230c get basis set handles; relativistic integral option only valid 1231c if both ibas and jbas are the ao basis. 1232c 1233 lbas = lc_bsh + BASIS_HANDLE_OFFSET 1234 sbas = sc_bsh + BASIS_HANDLE_OFFSET 1235 abas = ao_bsh + BASIS_HANDLE_OFFSET 1236 orel = ibas .eq. abas .and. jbas .eq. abas 1237 end if 1238c 1239c i shell 1240c 1241 ucont = (sf_ibs_cn2ucn(ish,ibas)) 1242c 1243c check for relativistic shell 1244c 1245 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then 1246 oirel = .true. 1247 isbas = sbas 1248 uconts = ao_to_ls(ucont) 1249 if (uconts .eq. 0) call errquit ( 1250 & 'int_1eefc1: no relativistic pointer',911, INT_ERR) 1251 if (nesc_1e_approx) then 1252 ibas = lbas 1253 ucont = uconts 1254 end if 1255 else 1256 uconts = ucont 1257 isbas = ibas 1258 end if 1259c 1260 Li = infbs_cont(CONT_TYPE ,ucont,ibas) 1261 i_prim = infbs_cont(CONT_NPRIM,ucont,ibas) 1262 i_gen = infbs_cont(CONT_NGEN ,ucont,ibas) 1263 i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 1264 i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas) 1265 i_cent = (sf_ibs_cn2ce(ish,ibas)) 1266 i_geom = ibs_geom(ibas) 1267 i_icfpS = infbs_cont(CONT_ICFP ,uconts,isbas) 1268c 1269c j shell 1270c 1271 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 1272c 1273c check for relativistic shell 1274c 1275 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then 1276 ojrel = .true. 1277 jsbas = sbas 1278 uconts = ao_to_ls(ucont) 1279 if (uconts .eq. 0) call errquit ( 1280 & 'int_1eefc1: no relativistic pointer',911, INT_ERR) 1281 if (nesc_1e_approx) then 1282 jbas = lbas 1283 ucont = uconts 1284 end if 1285 else 1286 uconts = ucont 1287 jsbas = jbas 1288 end if 1289c 1290 Lj = infbs_cont(CONT_TYPE ,ucont,jbas) 1291 j_prim = infbs_cont(CONT_NPRIM,ucont,jbas) 1292 j_gen = infbs_cont(CONT_NGEN ,ucont,jbas) 1293 j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas) 1294 j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas) 1295 j_cent = (sf_ibs_cn2ce(jsh,jbas)) 1296 j_geom = ibs_geom(jbas) 1297 j_icfpS = infbs_cont(CONT_ICFP ,uconts,jsbas) 1298c 1299 oNR = .not.(oirel.and.ojrel) 1300 orel = oirel.or.ojrel 1301c 1302 if (i_geom.ne.j_geom.and.WarnP.eq.0) then 1303 write(luout,*) 1304 & 'int_1eefc1: WARNING: possible geometry inconsistency' 1305 write(luout,*)'i_basis geometry handle:',i_geom 1306 write(luout,*)'j_basis geometry handle:',j_geom 1307 WarnP = 1 1308 endif 1309c 1310c ----- current cosmo restriction ... ----- 1311c 1312 if(i_geom.ne.j_geom) then 1313 write(luout,*) 'i_geom andr j_geom is noteq' 1314 call errquit('int_1eefc1: i_geom and j_geom must be the same', 1315 & 0,INT_ERR) 1316 endif 1317 1318 ohnd_ok = cando_hnd_1e(i_basis,ish,0) 1319 & .and. cando_hnd_1e(j_basis,jsh,0) 1320 & .and. (.not.geom_any_finuc (i_geom)) 1321 & .and. (.not.geom_any_finuc (j_geom)) 1322 onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0) 1323 onw1e_ok = cando_nw_1e(i_basis,ish,0) 1324 & .and. cando_nw_1e(j_basis,jsh,0) 1325c 1326c get external charges here (MV) 1327c ------------------------------ 1328 if(.not.geom_extbq_on()) 1329 > call errquit('int_1eefc1:no active bqs',0,0) 1330 bq_ncent = geom_extbq_ncenter() 1331 i_cbq = geom_extbq_coord() 1332 i_qbq = geom_extbq_charge() 1333 if(.not.ma_push_get(mt_dbl,bq_ncent,'xbq',h_xbq0,i_xbq0)) 1334 + call errquit( 'int_1eefc1',0,0) 1335 call dfill(bq_ncent,0.d0,dbl_mb(i_xbq0),1) 1336c 1337 if (orel) then 1338 call rel_onel ( 1339 & coords(1,i_cent,i_geom), 1340 & dbl_mb(mb_exndcf(i_iexp,ibas)), 1341 & dbl_mb(mb_exndcf(i_icfp,ibas)), 1342 & dbl_mb(mb_exndcf(i_icfpS,isbas)),i_prim,i_gen,Li, 1343 & coords(1,j_cent,j_geom), 1344 & dbl_mb(mb_exndcf(j_iexp,jbas)), 1345 & dbl_mb(mb_exndcf(j_icfp,jbas)), 1346 & dbl_mb(mb_exndcf(j_icfpS,jsbas)),j_prim,j_gen,Lj, 1347c & coords(1,1,cosmo_geom_efc),charge(1,cosmo_geom_efc), 1348c & geom_invnucexp(1,i_geom),ncenter(cosmo_geom_efc), 1349 & dbl_mb(i_cbq),dbl_mb(i_qbq), 1350 & dbl_mb(i_xbq0),bq_ncent, 1351 & S,T,V,lstv,doS,doT,doV,canAB,onw_ok,ohnd_ok,oNR,.false., 1352 & scr,lscr,rel_dbg,rel_typ) 1353 else if (ohnd_ok) then 1354 call hnd_stvint( 1355 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1356 & dbl_mb(mb_exndcf(i_icfp,ibas)), 1357 & i_prim, i_gen, Li, 1358 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1359 & dbl_mb(mb_exndcf(j_icfp,jbas)), 1360 & j_prim, j_gen, Lj, 1361c & coords(1,1,cosmo_geom_efc),charge(1,cosmo_geom_efc), 1362c & ncenter(cosmo_geom_efc), 1363 & dbl_mb(i_cbq),dbl_mb(i_qbq), 1364 & bq_ncent, 1365 & S,T,V,lstv,doS,doT,doV,scr,lscr) 1366c 1367 elseif (onw1e_ok) then 1368 call int_hf1sp( 1369 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1370 & dbl_mb(mb_exndcf(i_icfp,ibas)), 1371 & i_prim, i_gen, Li, i_cent, 1372 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1373 & dbl_mb(mb_exndcf(j_icfp,jbas)), 1374 & j_prim, j_gen, Lj, j_cent, 1375c & coords(1,1,cosmo_geom_efc),charge(1,cosmo_geom_efc), 1376c & geom_invnucexp(1,i_geom),ncenter(cosmo_geom_efc), 1377 & dbl_mb(i_cbq),dbl_mb(i_qbq), 1378 & dbl_mb(i_xbq0),bq_ncent, 1379 & S,T,V,lstv,doS,doT,doV,canAB,.false., 1380 & scr,lscr,'int_1eefc') 1381 elseif (onw_ok) then 1382 call hf1( 1383 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1384 & dbl_mb(mb_exndcf(i_icfp,ibas)), i_prim, i_gen, Li, 1385 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1386 & dbl_mb(mb_exndcf(j_icfp,jbas)), j_prim, j_gen, Lj, 1387c & coords(1,1,cosmo_geom_efc),charge(1,cosmo_geom_efc), 1388c & geom_invnucexp(1,i_geom),ncenter(cosmo_geom_efc), 1389 & dbl_mb(i_cbq),dbl_mb(i_qbq), 1390 & dbl_mb(i_xbq0),bq_ncent, 1391 & S,T,V,lstv,doS,doT,doV,canAB,.false., 1392 & scr,lscr) 1393 else 1394 call errquit('int_1eefc1: could not do hnd, sp or nw integrals', 1395 & 0, INT_ERR) 1396 endif 1397c 1398 if(.not.ma_pop_stack(h_xbq0)) 1399 + call errquit( 'int_1eefc1',0,0) 1400c 1401* We now have the cartesian integral block(s) (jlo:jhi,ilo:ihi) 1402* 1403 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 1404 if (.not.any_spherical) return 1405c 1406c ... reset general contractions for sp shells to 1 since they are handled 1407c as a block of 4. Since int_nbf_* arrays are set to the appropriate size. 1408c 1409 if (li.eq.-1) i_gen = 1 1410 if (lj.eq.-1) j_gen = 1 1411c 1412 if (bas_spherical(ibas).and.bas_spherical(jbas)) then 1413*... transform both i and j integrals 1414 i_nbf_x = int_nbf_x(Li) 1415 i_nbf_s = int_nbf_s(Li) 1416 j_nbf_x = int_nbf_x(Lj) 1417 j_nbf_s = int_nbf_s(Lj) 1418c 1419 if (doS) call spcart_tran1e(S,scr, 1420 & j_nbf_x,i_nbf_x,Lj,j_gen, 1421 & j_nbf_s,i_nbf_s,Li,i_gen, 1422 & .false.) 1423 if (doT) call spcart_tran1e(T,scr, 1424 & j_nbf_x,i_nbf_x,Lj,j_gen, 1425 & j_nbf_s,i_nbf_s,Li,i_gen, 1426 & .false.) 1427 if (doV) call spcart_tran1e(V,scr, 1428 & j_nbf_x,i_nbf_x,Lj,j_gen, 1429 & j_nbf_s,i_nbf_s,Li,i_gen, 1430 & .false.) 1431 else if (bas_spherical(ibas)) then 1432*.. transform on i component 1433 i_nbf_x = int_nbf_x(Li) 1434 i_nbf_s = int_nbf_s(Li) 1435 j_nbf_x = int_nbf_x(Lj) 1436 j_nbf_s = j_nbf_x 1437 if (doS) call spcart_tran1e(S,scr, 1438 & j_nbf_x,i_nbf_x,0,j_gen, 1439 & j_nbf_s,i_nbf_s,Li,i_gen, 1440 & .false.) 1441 if (doT) call spcart_tran1e(T,scr, 1442 & j_nbf_x,i_nbf_x,0,j_gen, 1443 & j_nbf_s,i_nbf_s,Li,i_gen, 1444 & .false.) 1445 if (doV) call spcart_tran1e(V,scr, 1446 & j_nbf_x,i_nbf_x,0,j_gen, 1447 & j_nbf_s,i_nbf_s,Li,i_gen, 1448 & .false.) 1449 else if (bas_spherical(jbas)) then 1450*.. transform on j component 1451 i_nbf_x = int_nbf_x(Li) 1452 i_nbf_s = i_nbf_x 1453 j_nbf_x = int_nbf_x(Lj) 1454 j_nbf_s = int_nbf_s(Lj) 1455 if (doS) call spcart_tran1e(S,scr, 1456 & j_nbf_x,i_nbf_x,Lj,j_gen, 1457 & j_nbf_s,i_nbf_s,0,i_gen, 1458 & .false.) 1459 if (doT) call spcart_tran1e(T,scr, 1460 & j_nbf_x,i_nbf_x,Lj,j_gen, 1461 & j_nbf_s,i_nbf_s,0,i_gen, 1462 & .false.) 1463 if (doV) call spcart_tran1e(V,scr, 1464 & j_nbf_x,i_nbf_x,Lj,j_gen, 1465 & j_nbf_s,i_nbf_s,0,i_gen, 1466 & .false.) 1467 else 1468 call errquit( 1469 & 'int_1eefc1: should never reach transform blocked else',911, 1470 & INT_ERR) 1471 endif 1472 return 1473c 1474 end 1475C> 1476C> @} 1477cc AJL/BEGIN/FDE 1478c---------------------------------------------------------------------- 1479C> \brief Compute the 1-electron Bq embedding potential 1480C> 1481C> This routine computes the 1 electron integrals \f$V\f$: 1482C> \f{eqnarray*}{ 1483C> V & = & 1484C({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\ 1485C> & = & 1486C\int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}} 1487C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 1488C> \f} 1489C> where \f$ q_\alpha \f$ is a point charge situated at position \f$ 1490Cr_\alpha \f$. 1491C> The Bq embedding point charges are accessed in int_1eefc1. 1492C> 1493 subroutine int_1efde1(i_basis,ish,j_basis,jsh,fde_basis,lscr,scr, 1494 & lpe,Pe) 1495#include "util.fh" 1496c 1497 integer i_basis !< [Input] basis set handle for ish 1498 integer ish !< [Input] i shell/contraction 1499 integer j_basis !< [Input] basis set handle for jsh 1500 integer jsh !< [Input] j shell/contraction 1501 integer fde_basis !< [Input] FDE basis set handle 1502 integer lscr !< [Input] length of scratch array 1503 double precision scr(lscr) !< [Scratch] scratch array 1504 integer lpe !< [Input] length of Pe buffer 1505 double precision Pe(lpe) !< [Output] potential energy integrals 1506c 1507c call int_1epe(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr, 1508c & scr,Pe,.false.,.false.,.true.,'pot') 1509 call int_1efde(i_basis,ish,j_basis,jsh,fde_basis,lscr,scr,lpe,scr, 1510 & scr,Pe,.false.,.false.,.true.) 1511c doS doT doV 1512 return 1513 end 1514c---------------------------------------------------------------------- 1515C> \brief Compute the 1-electron integrals but the potential is due to 1516C> Bq embedding charges 1517C> 1518C> This routine computes the 1 electron integrals \f$V\f$: 1519C> \f{eqnarray*}{ 1520C> V & = & 1521C({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\ 1522C> & = & 1523C\int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}} 1524C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1} 1525C> \f} 1526C> where \f$ q_\alpha \f$ is a point charge situated at position \f$ 1527Cr_\alpha \f$. 1528C> The location of the Bq embedding point charges is specified in a 1529C> Bq instance that has been activated elsewhere. 1530C> 1531 subroutine int_1efde(i_basis,ish,j_basis,jsh,fde_basis,lscr,scr, 1532 & lstv,S,T,V,doS,doT,doV) 1533 implicit none 1534#include "nwc_const.fh" 1535#include "errquit.fh" 1536#include "basP.fh" 1537#include "basdeclsP.fh" 1538#include "geomP.fh" 1539#include "geom.fh" 1540#include "geobasmapP.fh" 1541#include "mafdecls.fh" 1542#include "bas_exndcf_dec.fh" 1543#include "bas_ibs_dec.fh" 1544#include "int_nbf.fh" 1545#include "stdio.fh" 1546#include "apiP.fh" 1547#include "rel_nwc.fh" 1548#include "util.fh" 1549C 1550C The following include block contains the geometry handle for the 1551C geometry and charges of the solvent for COSMO. This is a clear 1552C violation 1553C of the integral API, but for now the cleanest solution 1554C 1555#include "prop.fh" 1556c::external subroutines used 1557c... errquit 1558c::functions 1559 logical cando_hnd_1e 1560 logical cando_nw_1e 1561 logical cando_nw 1562 logical int_chk_init 1563 logical int_chk_sh 1564 external int_chk_init 1565 external int_chk_sh 1566 external cando_hnd_1e 1567 external cando_nw_1e 1568 external cando_nw 1569c::passed 1570c:tex-\begin{verbatim} 1571 integer i_basis !< [Input] basis set handle for ish 1572 integer ish !< [Input] i shell/contraction 1573 integer j_basis !< [Input] basis set handle for jsh 1574 integer jsh !< [Input] j shell/contraction 1575 integer fde_basis !< [Input] FDE basis set handle 1576 integer lscr !< [Input] length of scratch array 1577 integer lstv !< [Input] length of integral buffer 1578 double precision scr(lscr) !< [Scratch] scratch array 1579 double precision S(lstv) !< [Output] overlap integrals 1580 double precision T(lstv) !< [Output] kinetic energy integrals 1581 double precision V(lstv) !< [Output] potential energy integrals 1582 logical doS !< [Input] flag for overlap integrals 1583 logical doT !< [Input] flag for kinetic energy integrals 1584 logical doV !< [Input] flag for potential energy integrals 1585c:tex-\end{verbatim} 1586c::local 1587 logical ohnd_ok, onw_ok, onw1e_ok 1588 logical shells_ok, orel, oirel, ojrel, oNR, canAB 1589 integer i_geom, j_geom, ibas, jbas, ucont, uconts 1590 integer lbas, sbas, abas, isbas, jsbas 1591 integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_icfpS 1592 integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_icfpS 1593c 1594 integer fde_bas, fde_geom 1595c 1596 logical any_spherical 1597 integer i_nbf_x, j_nbf_x 1598 integer i_nbf_s, j_nbf_s 1599c 1600 integer WarnP 1601 save WarnP 1602 data WarnP /0/ 1603c 1604#include "bas_exndcf_sfn.fh" 1605#include "bas_ibs_sfn.fh" 1606c 1607c check initialization and shells 1608c 1609 if (.not.int_chk_init('int_1efde')) 1610 & call errquit('int_1efde: int_init was not called' ,0, 1611 & INT_ERR) 1612c 1613 shells_ok = int_chk_sh(i_basis,ish) 1614 shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh) 1615 if (.not.shells_ok) 1616 & call errquit('int_1efde: invalid contraction/shell',0, 1617 & BASIS_ERR) 1618c 1619 ibas = i_basis + BASIS_HANDLE_OFFSET 1620 jbas = j_basis + BASIS_HANDLE_OFFSET 1621 fde_bas = fde_basis + BASIS_HANDLE_OFFSET 1622 oNR = .true. 1623 oirel = .false. 1624 ojrel = .false. 1625 orel = .false. 1626 canAB = .false. 1627c 1628 if (dyall_mod_dir) then 1629c 1630c get basis set handles; relativistic integral option only valid 1631c if both ibas and jbas are the ao basis. 1632c 1633 lbas = lc_bsh + BASIS_HANDLE_OFFSET 1634 sbas = sc_bsh + BASIS_HANDLE_OFFSET 1635 abas = ao_bsh + BASIS_HANDLE_OFFSET 1636 orel = ibas .eq. abas .and. jbas .eq. abas 1637 end if 1638c 1639c i shell 1640c 1641 ucont = (sf_ibs_cn2ucn(ish,ibas)) 1642c 1643c check for relativistic shell 1644c 1645 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then 1646 oirel = .true. 1647 isbas = sbas 1648 uconts = ao_to_ls(ucont) 1649 if (uconts .eq. 0) call errquit ( 1650 & 'int_1efde: no relativistic pointer',911, INT_ERR) 1651 if (nesc_1e_approx) then 1652 ibas = lbas 1653 ucont = uconts 1654 end if 1655 else 1656 uconts = ucont 1657 isbas = ibas 1658 end if 1659c 1660 Li = infbs_cont(CONT_TYPE ,ucont,ibas) 1661 i_prim = infbs_cont(CONT_NPRIM,ucont,ibas) 1662 i_gen = infbs_cont(CONT_NGEN ,ucont,ibas) 1663 i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 1664 i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas) 1665 i_cent = (sf_ibs_cn2ce(ish,ibas)) 1666 i_geom = ibs_geom(ibas) 1667 i_icfpS = infbs_cont(CONT_ICFP ,uconts,isbas) 1668c 1669c j shell 1670c 1671 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 1672c 1673c check for relativistic shell 1674c 1675 if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then 1676 ojrel = .true. 1677 jsbas = sbas 1678 uconts = ao_to_ls(ucont) 1679 if (uconts .eq. 0) call errquit ( 1680 & 'int_1efde: no relativistic pointer',911, INT_ERR) 1681 if (nesc_1e_approx) then 1682 jbas = lbas 1683 ucont = uconts 1684 end if 1685 else 1686 uconts = ucont 1687 jsbas = jbas 1688 end if 1689c 1690 Lj = infbs_cont(CONT_TYPE ,ucont,jbas) 1691 j_prim = infbs_cont(CONT_NPRIM,ucont,jbas) 1692 j_gen = infbs_cont(CONT_NGEN ,ucont,jbas) 1693 j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas) 1694 j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas) 1695 j_cent = (sf_ibs_cn2ce(jsh,jbas)) 1696 j_geom = ibs_geom(jbas) 1697 j_icfpS = infbs_cont(CONT_ICFP ,uconts,jsbas) 1698c 1699 oNR = .not.(oirel.and.ojrel) 1700 orel = oirel.or.ojrel 1701c 1702 if (i_geom.ne.j_geom.and.WarnP.eq.0) then 1703 write(luout,*) 1704 & 'int_1efde: WARNING: possible geometry inconsistency' 1705 write(luout,*)'i_basis geometry handle:',i_geom 1706 write(luout,*)'j_basis geometry handle:',j_geom 1707 WarnP = 1 1708 endif 1709c 1710c ----- current cosmo restriction ... ----- 1711c 1712 if(i_geom.ne.j_geom) then 1713 write(luout,*) 'i_geom and j_geom is not eq' 1714 call errquit('int_1efde: i_geom and j_geom must be the same', 1715 & 0,INT_ERR) 1716 endif 1717 1718 ohnd_ok = cando_hnd_1e(i_basis,ish,0) 1719 & .and. cando_hnd_1e(j_basis,jsh,0) 1720 & .and. (.not.geom_any_finuc (i_geom)) 1721 & .and. (.not.geom_any_finuc (j_geom)) 1722 onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0) 1723 onw1e_ok = cando_nw_1e(i_basis,ish,0) 1724 & .and. cando_nw_1e(j_basis,jsh,0) 1725c 1726c get external charges here (MV) 1727c ------------------------------ 1728c if(.not.geom_extbq_on()) 1729c > call errquit('int_1eefc1:no active bqs',0,0) 1730c bq_ncent = geom_extbq_ncenter() 1731c i_cbq = geom_extbq_coord() 1732c i_qbq = geom_extbq_charge() 1733c if(.not.ma_push_get(mt_dbl,bq_ncent,'xbq',h_xbq0,i_xbq0)) 1734c + call errquit( 'int_1eefc1',0,0) 1735c call dfill(bq_ncent,0.d0,dbl_mb(i_xbq0),1) 1736c 1737 fde_geom = ibs_geom(fde_bas) 1738 1739c write(*,*) 'FDE_GEOM ', fde_geom 1740c write(*,*) 'COORDS ', coords(1,1,fde_geom) 1741c write(*,*) 'CHARGE ', charge(1,fde_geom) 1742c write(*,*) 'INVNUCLEARCHARGE', geom_invnucexp(1,fde_geom) 1743c write(*,*) 'NCENTER ', ncenter(fde_geom) 1744 1745c call util_flush(luout) 1746 1747 if (orel) then 1748 call rel_onel ( 1749 & coords(1,i_cent,i_geom), 1750 & dbl_mb(mb_exndcf(i_iexp,ibas)), 1751 & dbl_mb(mb_exndcf(i_icfp,ibas)), 1752 & dbl_mb(mb_exndcf(i_icfpS,isbas)),i_prim,i_gen,Li, 1753 & coords(1,j_cent,j_geom), 1754 & dbl_mb(mb_exndcf(j_iexp,jbas)), 1755 & dbl_mb(mb_exndcf(j_icfp,jbas)), 1756 & dbl_mb(mb_exndcf(j_icfpS,jsbas)),j_prim,j_gen,Lj, 1757 & coords(1,1,fde_geom),charge(1,fde_geom), 1758 & geom_invnucexp(1,fde_geom),ncenter(fde_geom), 1759 & S,T,V,lstv,doS,doT,doV,canAB,onw_ok,ohnd_ok,oNR,.false., 1760 & scr,lscr,rel_dbg,rel_typ) 1761 else if (ohnd_ok) then 1762 call hnd_stvint( 1763 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1764 & dbl_mb(mb_exndcf(i_icfp,ibas)), 1765 & i_prim, i_gen, Li, 1766 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1767 & dbl_mb(mb_exndcf(j_icfp,jbas)), 1768 & j_prim, j_gen, Lj, 1769 & coords(1,1,fde_geom),charge(1,fde_geom), 1770 & ncenter(fde_geom), 1771 & S,T,V,lstv,doS,doT,doV,scr,lscr) 1772c 1773 elseif (onw1e_ok) then 1774 call int_hf1sp( 1775 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1776 & dbl_mb(mb_exndcf(i_icfp,ibas)), 1777 & i_prim, i_gen, Li, i_cent, 1778 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1779 & dbl_mb(mb_exndcf(j_icfp,jbas)), 1780 & j_prim, j_gen, Lj, j_cent, 1781 & coords(1,1,fde_geom),charge(1,fde_geom), 1782 & geom_invnucexp(1,fde_geom),ncenter(fde_geom), 1783 & S,T,V,lstv,doS,doT,doV,canAB,.false., 1784 & scr,lscr,'int_1eefc') 1785 elseif (onw_ok) then 1786 call hf1( 1787 & coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)), 1788 & dbl_mb(mb_exndcf(i_icfp,ibas)), i_prim, i_gen, Li, 1789 & coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)), 1790 & dbl_mb(mb_exndcf(j_icfp,jbas)), j_prim, j_gen, Lj, 1791 & coords(1,1,fde_geom),charge(1,fde_geom), 1792 & geom_invnucexp(1,fde_geom),ncenter(fde_geom), 1793 & S,T,V,lstv,doS,doT,doV,canAB,.false., 1794 & scr,lscr) 1795 else 1796 call errquit('int_1efde: could not do hnd, sp or nw integrals', 1797 & 0, INT_ERR) 1798 endif 1799c 1800c if(.not.ma_pop_stack(h_xbq0)) 1801c + call errquit( 'int_1eefc1',0,0) 1802c 1803* We now have the cartesian integral block(s) (jlo:jhi,ilo:ihi) 1804* 1805 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 1806 if (.not.any_spherical) return 1807c 1808c ... reset general contractions for sp shells to 1 since they are handled 1809c as a block of 4. Since int_nbf_* arrays are set to the appropriate size. 1810c 1811 if (li.eq.-1) i_gen = 1 1812 if (lj.eq.-1) j_gen = 1 1813c 1814 if (bas_spherical(ibas).and.bas_spherical(jbas)) then 1815*... transform both i and j integrals 1816 i_nbf_x = int_nbf_x(Li) 1817 i_nbf_s = int_nbf_s(Li) 1818 j_nbf_x = int_nbf_x(Lj) 1819 j_nbf_s = int_nbf_s(Lj) 1820c 1821 if (doS) call spcart_tran1e(S,scr, 1822 & j_nbf_x,i_nbf_x,Lj,j_gen, 1823 & j_nbf_s,i_nbf_s,Li,i_gen, 1824 & .false.) 1825 if (doT) call spcart_tran1e(T,scr, 1826 & j_nbf_x,i_nbf_x,Lj,j_gen, 1827 & j_nbf_s,i_nbf_s,Li,i_gen, 1828 & .false.) 1829 if (doV) call spcart_tran1e(V,scr, 1830 & j_nbf_x,i_nbf_x,Lj,j_gen, 1831 & j_nbf_s,i_nbf_s,Li,i_gen, 1832 & .false.) 1833 else if (bas_spherical(ibas)) then 1834*.. transform on i component 1835 i_nbf_x = int_nbf_x(Li) 1836 i_nbf_s = int_nbf_s(Li) 1837 j_nbf_x = int_nbf_x(Lj) 1838 j_nbf_s = j_nbf_x 1839 if (doS) call spcart_tran1e(S,scr, 1840 & j_nbf_x,i_nbf_x,0,j_gen, 1841 & j_nbf_s,i_nbf_s,Li,i_gen, 1842 & .false.) 1843 if (doT) call spcart_tran1e(T,scr, 1844 & j_nbf_x,i_nbf_x,0,j_gen, 1845 & j_nbf_s,i_nbf_s,Li,i_gen, 1846 & .false.) 1847 if (doV) call spcart_tran1e(V,scr, 1848 & j_nbf_x,i_nbf_x,0,j_gen, 1849 & j_nbf_s,i_nbf_s,Li,i_gen, 1850 & .false.) 1851 else if (bas_spherical(jbas)) then 1852*.. transform on j component 1853 i_nbf_x = int_nbf_x(Li) 1854 i_nbf_s = i_nbf_x 1855 j_nbf_x = int_nbf_x(Lj) 1856 j_nbf_s = int_nbf_s(Lj) 1857 if (doS) call spcart_tran1e(S,scr, 1858 & j_nbf_x,i_nbf_x,Lj,j_gen, 1859 & j_nbf_s,i_nbf_s,0,i_gen, 1860 & .false.) 1861 if (doT) call spcart_tran1e(T,scr, 1862 & j_nbf_x,i_nbf_x,Lj,j_gen, 1863 & j_nbf_s,i_nbf_s,0,i_gen, 1864 & .false.) 1865 if (doV) call spcart_tran1e(V,scr, 1866 & j_nbf_x,i_nbf_x,Lj,j_gen, 1867 & j_nbf_s,i_nbf_s,0,i_gen, 1868 & .false.) 1869 else 1870 call errquit( 1871 & 'int_1efde: should never reach transform blocked else',911, 1872 & INT_ERR) 1873 endif 1874 return 1875c 1876 end 1877cc AJL/End 1878C> 1879C> @} 1880