1C> \ingroup nwint 2C> @{ 3C> 4C> \brief Compute the multipole integral derivatives 5C> 6C> Compute the multipole integral derivatives of the general form 7C> \f{eqnarray*}{ 8C> \frac{\partial (\mu|M_A|\nu)}{\partial R_x} &=& 9C> \int \frac{\partial [g_\mu(X_\mu,r_1)(R_A-r_1)_x^{n_x}(R_A-r_1)_y^{n_y}(R_A-r_1)_z^{n_z}g_\nu(X_\nu,r_1)]}{\partial X_x}dr_1 10C> \f} 11C> where the output buffer is logically organized as 12C> `MP(mrange,jlo:jhi,ilo:ihi,ncoord,natom)` where 13C> 14C> * `mrange` is the range of multipoles 1:((l+1)*(l+2)/2) 15C> 16C> * `jlo:jhi` is the range of basis functions associated with `jsh` 17C> 18C> * `ilo:ihi` is the range of basis functions associated with `ish` 19C> 20C> * `ncoord` represents the 3 Cartesian coordinates of the derivative 21C> 22C> * `natom` represents the 3 "atoms" involved in the integral, where 23C> 1 is the center associated with `ish`, 2 is the center associate with 24C> `jsh`, and 3 is associated with the center of the multipole. 25C> 26C> Integral derivatives are returned in shell blocks of <L|ish|jsh> L=lval 27C> one block for the given L value. EACH block repeated for 9 (xyz,atoms) 28C> for ish = d and Lval = 1 and jsh = p you would get: 29C> 30C> * (6*3*3)*3*3=486 integral derivatives 31C> 32C> order would be 33C> \code 34C> <xx|x|x> <xx|y|x> <xx|z|x> ( 1- 3)(x,atom1) 35C> <xx|x|y> <xx|y|y> <xx|z|y> ( 4- 6)(x,atom1) 36C> <xx|x|z> <xx|y|z> <xx|z|z> ( 7- 9)(x,atom1) 37C> <xy|x|x> <xy|y|x> <xy|z|x> (10-12)(x,atom1) 38C> <xy|x|y> <xy|y|y> <xy|z|y> (13-15)(x,atom1) 39C> <xy|x|z> <xy|y|z> <xy|z|z> (16-18)(x,atom1) 40C> <xz|x|x> <xz|y|x> <xz|z|x> (19-21)(x,atom1) 41C> <xz|x|y> <xz|y|y> <xz|z|y> (22-24)(x,atom1) 42C> <xz|x|z> <xz|y|z> <xz|z|z> (25-27)(x,atom1) 43C> <yy|x|x> <yy|y|x> <yy|z|x> (28-30)(x,atom1) 44C> <yy|x|y> <yy|y|y> <yy|z|y> (31-33)(x,atom1) 45C> <yy|x|z> <yy|y|z> <yy|z|z> (34-36)(x,atom1) 46C> <yz|x|x> <yz|y|x> <yz|z|x> (37-39)(x,atom1) 47C> <yz|x|y> <yz|y|y> <yz|z|y> (40-42)(x,atom1) 48C> <yz|x|z> <yz|y|z> <yz|z|z> (43-45)(x,atom1) 49C> <zz|x|x> <zz|y|x> <zz|z|x> (46-48)(x,atom1) 50C> <zz|x|y> <zz|y|y> <zz|z|y> (49-51)(x,atom1) 51C> <zz|x|z> <zz|y|z> <zz|z|z> (52-54)(x,atom1) 52C> repeat above for (y,atom1), (z,atom1) 53C> repeat above for atom2 and multipole center 54C> \endcode 55C> 56C> 57C> 58C> For ish = p and Lval = 2 and jsh = p you would get: 59C> 60C> * (3*6*3)*3*3 = 486 integral derivatives 61C> 62C> The order would be 63C> \code 64C> <x|xx|x> <x|xy|x> <x|xz|x> <x|yy|x> <x|yz|x> <x|zz|x> ( 1- 6)(x,atom1) 65C> <x|xx|y> <x|xy|y> <x|xz|y> <x|yy|y> <x|yz|y> <x|zz|y> ( 7-12)(x,atom1) 66C> <x|xx|z> <x|xy|z> <x|xz|z> <x|yy|z> <x|yz|z> <x|zz|z> (13-18)(x,atom1) 67C> <y|xx|x> <y|xy|x> <y|xz|x> <y|yy|x> <y|yz|x> <y|zz|x> (19-24)(x,atom1) 68C> <y|xx|y> <y|xy|y> <y|xz|y> <y|yy|y> <y|yz|y> <y|zz|y> (25-30)(x,atom1) 69C> <y|xx|z> <y|xy|z> <y|xz|z> <y|yy|z> <y|yz|z> <y|zz|z> (31-36)(x,atom1) 70C> <z|xx|x> <z|xy|x> <z|xz|x> <z|yy|x> <z|yz|x> <z|zz|x> (37-42)(x,atom1) 71C> <z|xx|y> <z|xy|y> <z|xz|y> <z|yy|y> <z|yz|y> <z|zz|y> (43-48)(x,atom1) 72C> <z|xx|z> <z|xy|z> <z|xz|z> <z|yy|z> <z|yz|z> <z|zz|z> (49-54)(x,atom1) 73C> repeat above for (y,atom1), (z,atom1) 74C> repeat above for atom2 and multipole center 75C> \endcode 76C> 77C> For ish = s and lval = 4 and jsh = p you would get: 78C> 79C> * (1*15*3)*3*3 = 405 integral derivatives 80C> 81C> and the order would be 82C> \code 83C> <s|xxxx|x> <s|xxxy|x> <s|xxxz|x> <s|xxyy|x> <s|xxyz|x> <s|xxzz|x> ( 1- 6)(x,atom1) 84C> <s|xyyy|x> <s|xyyz|x> <s|xyzz|x> <s|xzzz|x> <s|yyyy|x> <s|yyyz|x> ( 7-12)(x,atom1) 85C> <s|yyzz|x> <s|yzzz|x> <s|zzzz|x> (13-15)(x,atom1) 86C> <s|xxxx|y> <s|xxxy|y> <s|xxxz|y> <s|xxyy|y> <s|xxyz|y> <s|xxzz|y> (16-21)(x,atom1) 87C> <s|xyyy|y> <s|xyyz|y> <s|xyzz|y> <s|xzzz|y> <s|yyyy|y> <s|yyyz|y> (22-27)(x,atom1) 88C> <s|yyzz|y> <s|yzzz|y> <s|zzzz|y> (28-30)(x,atom1) 89C> <s|xxxx|z> <s|xxxy|z> <s|xxxz|z> <s|xxyy|z> <s|xxyz|z> <s|xxzz|z> (31-36)(x,atom1) 90C> <s|xyyy|z> <s|xyyz|z> <s|xyzz|z> <s|xzzz|z> <s|yyyy|z> <s|yyyz|z> (37-42)(x,atom1) 91C> <s|yyzz|z> <s|yzzz|z> <s|zzzz|z> (43-45)(x,atom1) 92C> repeat above for (y,atom1), (z,atom1) 93C> repeat above for atom2 and multipole center 94C> \endcode 95C> 96 subroutine intd_mpolel(i_basis, ish, j_basis, jsh, 97 & lval, centerl, 98 & lscr, scr, lmpint, MP, num_mpint, 99 & idatom) 100* 101* $Id$ 102* 103c 104c routine to compute multipole integral derivatives at a given lvalue 105c 106c The general form is <shell|pole|shell|3|3> 107c 108c the returned buffer is logically: 109c (mpole range, jlo:jhi, ilo:ihi, 3(xyz), 3(atom)) 110c where mpole range is 1:((lval+1)*(lval+2)/2) 111c 3(xyz) is the x,y,z coordinate derivative for the atom index 112c 3(atom) is at most three centers one is where ish exists 113c two is where jsh exists 114c three is the center of the multipole 115c 116c Integral derivatives are returned in shell blocks of <L|ish|jsh> L=lval 117c one block for the given L value. EACH block repeated for by 9 (xyz,atoms) 118c for ish = d and Lval = 1 and jsh = p you would get: 119c (6*3*3)*3*3=486 integral derivatives 120c order would be 121c <xx|x|x> <xx|y|x> <xx|z|x> ( 1- 3)(x,atom1) 122c <xx|x|y> <xx|y|y> <xx|z|y> ( 4- 6)(x,atom1) 123c <xx|x|z> <xx|y|z> <xx|z|z> ( 7- 9)(x,atom1) 124c <xy|x|x> <xy|y|x> <xy|z|x> (10-12)(x,atom1) 125c <xy|x|y> <xy|y|y> <xy|z|y> (13-15)(x,atom1) 126c <xy|x|z> <xy|y|z> <xy|z|z> (16-18)(x,atom1) 127c <xz|x|x> <xz|y|x> <xz|z|x> (19-21)(x,atom1) 128c <xz|x|y> <xz|y|y> <xz|z|y> (22-24)(x,atom1) 129c <xz|x|z> <xz|y|z> <xz|z|z> (25-27)(x,atom1) 130c <yy|x|x> <yy|y|x> <yy|z|x> (28-30)(x,atom1) 131c <yy|x|y> <yy|y|y> <yy|z|y> (31-33)(x,atom1) 132c <yy|x|z> <yy|y|z> <yy|z|z> (34-36)(x,atom1) 133c <yz|x|x> <yz|y|x> <yz|z|x> (37-39)(x,atom1) 134c <yz|x|y> <yz|y|y> <yz|z|y> (40-42)(x,atom1) 135c <yz|x|z> <yz|y|z> <yz|z|z> (43-45)(x,atom1) 136c <zz|x|x> <zz|y|x> <zz|z|x> (46-48)(x,atom1) 137c <zz|x|y> <zz|y|y> <zz|z|y> (49-51)(x,atom1) 138c <zz|x|z> <zz|y|z> <zz|z|z> (52-54)(x,atom1) 139c repeat above for (y,atom1), (z,atom1) 140c repeat above for atom2 and multipole center 141c 142c 143c 144c for ish = p and Lval = 2 and jsh = p you would get: 145c (3*6*3)*3*3 = 486 integral derivatives 146c order would be 147c <x|xx|x> <x|xy|x> <x|xz|x> <x|yy|x> <x|yz|x> <x|zz|x> ( 1- 6)(x,atom1) 148c <x|xx|y> <x|xy|y> <x|xz|y> <x|yy|y> <x|yz|y> <x|zz|y> ( 7-12)(x,atom1) 149c <x|xx|z> <x|xy|z> <x|xz|z> <x|yy|z> <x|yz|z> <x|zz|z> (13-18)(x,atom1) 150c <y|xx|x> <y|xy|x> <y|xz|x> <y|yy|x> <y|yz|x> <y|zz|x> (19-24)(x,atom1) 151c <y|xx|y> <y|xy|y> <y|xz|y> <y|yy|y> <y|yz|y> <y|zz|y> (25-30)(x,atom1) 152c <y|xx|z> <y|xy|z> <y|xz|z> <y|yy|z> <y|yz|z> <y|zz|z> (31-36)(x,atom1) 153c <z|xx|x> <z|xy|x> <z|xz|x> <z|yy|x> <z|yz|x> <z|zz|x> (37-42)(x,atom1) 154c <z|xx|y> <z|xy|y> <z|xz|y> <z|yy|y> <z|yz|y> <z|zz|y> (43-48)(x,atom1) 155c <z|xx|z> <z|xy|z> <z|xz|z> <z|yy|z> <z|yz|z> <z|zz|z> (49-54)(x,atom1) 156c repeat above for (y,atom1), (z,atom1) 157c repeat above for atom2 and multipole center 158c 159c for ish = s and lval = 4 and jsh = p you would get: 160c (1*15*3)*3*3 = 405 integral derivatives 161c <s|xxxx|x> <s|xxxy|x> <s|xxxz|x> <s|xxyy|x> <s|xxyz|x> <s|xxzz|x> ( 1- 6)(x,atom1) 162c <s|xyyy|x> <s|xyyz|x> <s|xyzz|x> <s|xzzz|x> <s|yyyy|x> <s|yyyz|x> ( 7-12)(x,atom1) 163c <s|yyzz|x> <s|yzzz|x> <s|zzzz|x> (13-15)(x,atom1) 164c <s|xxxx|y> <s|xxxy|y> <s|xxxz|y> <s|xxyy|y> <s|xxyz|y> <s|xxzz|y> (16-21)(x,atom1) 165c <s|xyyy|y> <s|xyyz|y> <s|xyzz|y> <s|xzzz|y> <s|yyyy|y> <s|yyyz|y> (22-27)(x,atom1) 166c <s|yyzz|y> <s|yzzz|y> <s|zzzz|y> (28-30)(x,atom1) 167c <s|xxxx|z> <s|xxxy|z> <s|xxxz|z> <s|xxyy|z> <s|xxyz|z> <s|xxzz|z> (31-36)(x,atom1) 168c <s|xyyy|z> <s|xyyz|z> <s|xyzz|z> <s|xzzz|z> <s|yyyy|z> <s|yyyz|z> (37-42)(x,atom1) 169c <s|yyzz|z> <s|yzzz|z> <s|zzzz|z> (43-45)(x,atom1) 170c repeat above for (y,atom1), (z,atom1) 171c repeat above for atom2 and multipole center 172c 173c 174 implicit none 175#include "apiP.fh" 176#include "errquit.fh" 177#include "nwc_const.fh" 178#include "basP.fh" 179#include "basdeclsP.fh" 180#include "geobasmapP.fh" 181#include "geomP.fh" 182#include "stdio.fh" 183#include "mafdecls.fh" 184#include "bas_exndcf_dec.fh" 185#include "bas_ibs_dec.fh" 186#include "int_nbf.fh" 187c 188c::functions 189 logical int_chk_init 190 integer int_nint_cart 191 external int_chk_init 192 external int_nint_cart 193c::passed 194 integer i_basis !< [Input] basis set handle for ish 195 integer ish !< [Input] i shell/contraction 196 integer j_basis !< [Input] basis set handle for jsh 197 integer jsh !< [Input] j shell/contraction 198 integer lval !< [Input] maximum lvalue for 199*.......................................... multipole integrals 200*.......................................... in this batch 201 double precision centerl(3) !< [Input] coordinates of multipole 202 integer lscr !< [Input] length of scratch array 203 double precision scr(lscr) !< [Input] scratch array 204 integer lmpint !< [Input] length of multipole 205*.......................................... integrals array 206 double precision MP(lmpint) !< [Output] multipole integrals 207 integer num_mpint !< [Output] number of multipole integrals 208 integer idatom(3) !< [Output] array identifying centers for derivatives 209c ! e.g., the first nint*3 derivatives go to center idatom(1) 210c ! the second nint*3 derivatives go to center idatom(2) 211c::local 212 logical shells_ok 213 integer ibas, Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 214 integer jbas, Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 215 integer ucont 216 integer l_int, ij_int, num_int, num_intd 217 integer lpole 218 integer ioffj 219 logical any_spherical 220 logical inline_chk_sh 221c 222 integer WarnP 223 save WarnP 224 data WarnP /0/ 225c 226#include "bas_exndcf_sfn.fh" 227#include "bas_ibs_sfn.fh" 228c 229c... statement function for int_chk_sh 230 inline_chk_sh(ibas,ish) = 231 $ ((ish.gt.0) .and. (ish.le.ncont_tot_gb(ibas))) 232c 233c check initialization 234c 235 if (.not.int_chk_init('intd_mpolel')) 236 & call errquit('intd_mpolel: int_init was not called' ,0, 237 & INT_ERR) 238c 239c check if gencon/sp shells 240c 241 call int_nogencont_check(i_basis,'intd_mpolel:i_basis') 242 call int_nogencont_check(j_basis,'intd_mpolel:j_basis') 243 call int_nospshell_check(i_basis,'intd_mpolel:i_basis') 244 call int_nospshell_check(j_basis,'intd_mpolel:j_basis') 245c 246 ibas = i_basis + BASIS_HANDLE_OFFSET 247 jbas = j_basis + BASIS_HANDLE_OFFSET 248c 249 shells_ok = inline_chk_sh(ibas,ish) 250 shells_ok = shells_ok .and. inline_chk_sh(jbas,jsh) 251 if (.not. shells_ok) 252 & call errquit('intd_mpolel: invalid contraction/shell',0, 253 & BASIS_ERR) 254c 255*** set defNxyz such that it can handle the maximum multi-pole 256c 257 lpole = lval/4 + 1 258 lpole = lpole + 1 ! for derivative 259 call defNxyz(lpole) 260c 261 ucont = (sf_ibs_cn2ucn(ish,ibas)) 262 Li = infbs_cont(CONT_TYPE ,ucont,ibas) 263 i_prim = infbs_cont(CONT_NPRIM,ucont,ibas) 264 i_gen = infbs_cont(CONT_NGEN ,ucont,ibas) 265 i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 266 i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas) 267 i_cent = (sf_ibs_cn2ce(ish,ibas)) 268 i_geom = ibs_geom(ibas) 269c 270 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 271 Lj = infbs_cont(CONT_TYPE ,ucont,jbas) 272 j_prim = infbs_cont(CONT_NPRIM,ucont,jbas) 273 j_gen = infbs_cont(CONT_NGEN ,ucont,jbas) 274 j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas) 275 j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas) 276 j_cent = (sf_ibs_cn2ce(jsh,jbas)) 277 j_geom = ibs_geom(jbas) 278c 279 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 280c 281 if (i_geom.ne.j_geom.and.WarnP.eq.0) then 282 write(luout,*)'intd_mpolel: WARNING: possible geometry', 283 & ' inconsistency' 284 write(luout,*)'i_basis geometry handle:',i_geom 285 write(luout,*)'j_basis geometry handle:',j_geom 286 WarnP = 1 287 endif 288c 289 if (i_gen.gt.1 .or. j_gen.gt.1) then 290 write(luout,*) 291 & ' hf3ois does not handle general contractions yet' 292 call errquit('intd_mpolel: general contraction error ',911, 293 & INT_ERR) 294 endif 295c 296 l_int = (lval+1)*(lval+2)/2 297 ij_int = int_nint_cart(i_basis, ish, j_basis, jsh, 0,0, 0,0) 298 num_int = l_int*ij_int 299 num_intd = num_int*3*3 300 if (num_intd.gt.lmpint) then 301 write(luout,*)' intd_mpolel: lmpint = ',lmpint 302 write(luout,*)' intd_mpolel: num_intd = ',num_intd 303 call errquit('intd_mpolel: lmpint too small ',911, INT_ERR) 304 endif 305 if ((lval.eq.0).and.(i_cent.eq.j_cent)) then 306 call dcopy(num_intd,0.0d00,0,MP,1) 307 call ifill(3,0,idatom,1) 308 return 309 endif 310c 311 call hfd3ois( 312 & coords(1,i_cent,i_geom), 313 & dbl_mb(mb_exndcf(i_iexp,ibas)), 314 & dbl_mb(mb_exndcf(i_icfp,ibas)), 315 & i_prim, 1, Li, 316 & coords(1,j_cent,j_geom), 317 & dbl_mb(mb_exndcf(j_iexp,jbas)), 318 & dbl_mb(mb_exndcf(j_icfp,jbas)), 319 & j_prim, 1, Lj, 320 & centerl, 321 & 0.0d00, 322 & 1.0d00, 323 & 1, 1, lval, 324c.....................DryRun 325 & MP,num_int,.false.,scr,lscr) 326 num_mpint = num_intd 327 if (i_cent.eq.j_cent) then 328 idatom(1) = i_cent 329 idatom(2) = 0 330 idatom(3) = -3 331 ioffj = num_int*3 + 1 332 call daxpy(num_int*3,1.0d00,MP(ioffj),1,MP,1) 333 else 334 idatom(1) = i_cent 335 idatom(2) = j_cent 336 idatom(3) = -3 337 endif 338 if (any_spherical) then 339 if (Li.eq.-1) i_gen = 1 340 if (Lj.eq.-1) j_gen = 1 341 call spcart_3cBtran(MP,scr,lscr, 342 & int_nbf_x(Li),int_nbf_s(Li),Li,i_gen,bas_spherical(ibas), 343 & int_nbf_x(Lj),int_nbf_s(Lj),Lj,j_gen,bas_spherical(jbas), 344 & int_nbf_x(lval),int_nbf_x(lval),lval,1,.false., 345 & 9,.false.) 346 endif 347 end 348C> @} 349