1C> \ingroup nwint 2C> @{ 3C> 4C> \brief Compute 1-electron Gaussian periodic multipole 5C> integral derivatives for \f$l\f$ from \f$0\f$ to \f$l_{max}\f$ 6C> 7C> See [1] for details. 8C> 9C> [1] JE Jaffe, AC Hess, 10C> <i>"Gaussian basis density functional theory for systems 11C> periodic in two or three dimensions: Energy and forces"</i>, 12C> J.Chem.Phys. <b>105</b>, 10983-10998 (1996), DOI: 13C> <a href="https://doi.org/10.1063/1.472866"> 14C> 10.1063/1.472866</a> 15C> 16 subroutine intpd_mpolel(i_basis, ish, j_basis, jsh, R, 17 & lval, centerl, 18 & lscr, scr, lmpint, MP, num_mpint, 19 & idatom) 20* 21* $Id$ 22* 23c 24c routine to compute multipole integral derivatives at a given lvalue 25c with the jsh translated by the fractional coordinate vector R. 26c The general form is <shell|pole|shell|3|3> 27c 28c the returned buffer is logically: 29c (mpole range, jlo:jhi, ilo:ihi, 3(xyz), 3(atom)) 30c where mpole range is 1:((lval+1)*(lval+2)/2) 31c 3(xyz) is the x,y,z coordinate derivative for the atom index 32c 3(atom) is at most three centers one is where ish exists 33c two is where jsh exists 34c three is the center of the multipole 35c 36c Integral derivatives are returned in shell blocks of <L|ish|jsh> L=lval 37c one block for the given L value. EACH block repeated for by 9 (xyz,atoms) 38c for ish = d and Lval = 1 and jsh = p you would get: 39c (6*3*3)*3*3=486 integral derivatives 40c order would be 41c <xx|x|x> <xx|y|x> <xx|z|x> ( 1- 3)(x,atom1) 42c <xx|x|y> <xx|y|y> <xx|z|y> ( 4- 6)(x,atom1) 43c <xx|x|z> <xx|y|z> <xx|z|z> ( 7- 9)(x,atom1) 44c <xy|x|x> <xy|y|x> <xy|z|x> (10-12)(x,atom1) 45c <xy|x|y> <xy|y|y> <xy|z|y> (13-15)(x,atom1) 46c <xy|x|z> <xy|y|z> <xy|z|z> (16-18)(x,atom1) 47c <xz|x|x> <xz|y|x> <xz|z|x> (19-21)(x,atom1) 48c <xz|x|y> <xz|y|y> <xz|z|y> (22-24)(x,atom1) 49c <xz|x|z> <xz|y|z> <xz|z|z> (25-27)(x,atom1) 50c <yy|x|x> <yy|y|x> <yy|z|x> (28-30)(x,atom1) 51c <yy|x|y> <yy|y|y> <yy|z|y> (31-33)(x,atom1) 52c <yy|x|z> <yy|y|z> <yy|z|z> (34-36)(x,atom1) 53c <yz|x|x> <yz|y|x> <yz|z|x> (37-39)(x,atom1) 54c <yz|x|y> <yz|y|y> <yz|z|y> (40-42)(x,atom1) 55c <yz|x|z> <yz|y|z> <yz|z|z> (43-45)(x,atom1) 56c <zz|x|x> <zz|y|x> <zz|z|x> (46-48)(x,atom1) 57c <zz|x|y> <zz|y|y> <zz|z|y> (49-51)(x,atom1) 58c <zz|x|z> <zz|y|z> <zz|z|z> (52-54)(x,atom1) 59c repeat above for (y,atom1), (z,atom1) 60c repeat above for atom2 and multipole center 61c 62c 63c 64c for ish = p and Lval = 2 and jsh = p you would get: 65c (3*6*3)*3*3 = 486 integral derivatives 66c order would be 67c <x|xx|x> <x|xy|x> <x|xz|x> <x|yy|x> <x|yz|x> <x|zz|x> ( 1- 6)(x,atom1) 68c <x|xx|y> <x|xy|y> <x|xz|y> <x|yy|y> <x|yz|y> <x|zz|y> ( 7-12)(x,atom1) 69c <x|xx|z> <x|xy|z> <x|xz|z> <x|yy|z> <x|yz|z> <x|zz|z> (13-18)(x,atom1) 70c <y|xx|x> <y|xy|x> <y|xz|x> <y|yy|x> <y|yz|x> <y|zz|x> (19-24)(x,atom1) 71c <y|xx|y> <y|xy|y> <y|xz|y> <y|yy|y> <y|yz|y> <y|zz|y> (25-30)(x,atom1) 72c <y|xx|z> <y|xy|z> <y|xz|z> <y|yy|z> <y|yz|z> <y|zz|z> (31-36)(x,atom1) 73c <z|xx|x> <z|xy|x> <z|xz|x> <z|yy|x> <z|yz|x> <z|zz|x> (37-42)(x,atom1) 74c <z|xx|y> <z|xy|y> <z|xz|y> <z|yy|y> <z|yz|y> <z|zz|y> (43-48)(x,atom1) 75c <z|xx|z> <z|xy|z> <z|xz|z> <z|yy|z> <z|yz|z> <z|zz|z> (49-54)(x,atom1) 76c repeat above for (y,atom1), (z,atom1) 77c repeat above for atom2 and multipole center 78c 79c for ish = s and lval = 4 and jsh = p you would get: 80c (1*15*3)*3*3 = 405 integral derivatives 81c <s|xxxx|x> <s|xxxy|x> <s|xxxz|x> <s|xxyy|x> <s|xxyz|x> <s|xxzz|x> ( 1- 6)(x,atom1) 82c <s|xyyy|x> <s|xyyz|x> <s|xyzz|x> <s|xzzz|x> <s|yyyy|x> <s|yyyz|x> ( 7-12)(x,atom1) 83c <s|yyzz|x> <s|yzzz|x> <s|zzzz|x> (13-15)(x,atom1) 84c <s|xxxx|y> <s|xxxy|y> <s|xxxz|y> <s|xxyy|y> <s|xxyz|y> <s|xxzz|y> (16-21)(x,atom1) 85c <s|xyyy|y> <s|xyyz|y> <s|xyzz|y> <s|xzzz|y> <s|yyyy|y> <s|yyyz|y> (22-27)(x,atom1) 86c <s|yyzz|y> <s|yzzz|y> <s|zzzz|y> (28-30)(x,atom1) 87c <s|xxxx|z> <s|xxxy|z> <s|xxxz|z> <s|xxyy|z> <s|xxyz|z> <s|xxzz|z> (31-36)(x,atom1) 88c <s|xyyy|z> <s|xyyz|z> <s|xyzz|z> <s|xzzz|z> <s|yyyy|z> <s|yyyz|z> (37-42)(x,atom1) 89c <s|yyzz|z> <s|yzzz|z> <s|zzzz|z> (43-45)(x,atom1) 90c repeat above for (y,atom1), (z,atom1) 91c repeat above for atom2 and multipole center 92c 93c 94 implicit none 95#include "apiP.fh" 96#include "errquit.fh" 97#include "nwc_const.fh" 98#include "basP.fh" 99#include "basdeclsP.fh" 100#include "geobasmapP.fh" 101#include "geomP.fh" 102#include "stdio.fh" 103#include "mafdecls.fh" 104#include "bas_exndcf_dec.fh" 105#include "bas_ibs_dec.fh" 106#include "int_nbf.fh" 107c 108c::functions 109 logical int_chk_init 110 integer int_nint_cart 111 external int_chk_init 112 external int_nint_cart 113c::passed 114 integer i_basis !< [Input] basis set handle for ish 115 integer ish !< [Input] i shell/contraction 116 integer j_basis !< [Input] basis set handle for jsh 117 integer jsh !< [Input] j shell/contraction 118c... translation vectors are in fractional coordinates 119 double precision R(3) !< [Input] translation vec on j cont. 120 integer lval !< [Input] maximum lvalue for 121 !< multipole integrals in this batch 122 double precision centerl(3) !< [Input] coordinates of multipole 123 integer lscr !< [Input] length of scratch array 124 double precision scr(lscr) !< [Input] scratch array 125 integer lmpint !< [Input] length of multipole 126 !< integrals array 127 double precision MP(lmpint) !< [Output] multipole integrals 128 integer num_mpint !< [Output] number of multipole integrals 129 integer idatom(3) !< [Output] array identifying centers for derivatives 130c ! e.g., the first nint*3 derivatives go to center idatom(1) 131c ! the second nint*3 derivatives go to center idatom(2) 132c::local 133 logical shells_ok 134 integer ibas, Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 135 integer jbas, Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 136 integer ucont 137 integer l_int, ij_int, num_int, num_intd 138 integer lpole 139 integer ioffj 140 double precision xyz_new(3) 141 logical any_spherical 142 logical inline_chk_sh 143c 144 integer WarnP 145 save WarnP 146 data WarnP /0/ 147c 148#include "bas_exndcf_sfn.fh" 149#include "bas_ibs_sfn.fh" 150c 151c... statement function for int_chk_sh 152 inline_chk_sh(ibas,ish) = 153 $ ((ish.gt.0) .and. (ish.le.ncont_tot_gb(ibas))) 154c 155c check initialization 156c 157 if (.not.int_chk_init('intpd_mpolel')) 158 & call errquit('intpd_mpolel: int_init was not called' ,0, 159 & INT_ERR) 160c 161c check if gencon/sp shells 162c 163 call int_nogencont_check(i_basis,'intpd_mpolel:i_basis') 164 call int_nogencont_check(j_basis,'intpd_mpolel:j_basis') 165 call int_nospshell_check(i_basis,'intpd_mpolel:i_basis') 166 call int_nospshell_check(j_basis,'intpd_mpolel:j_basis') 167c 168 ibas = i_basis + BASIS_HANDLE_OFFSET 169 jbas = j_basis + BASIS_HANDLE_OFFSET 170c 171 shells_ok = inline_chk_sh(ibas,ish) 172 shells_ok = shells_ok .and. inline_chk_sh(jbas,jsh) 173 if (.not. shells_ok) 174 & call errquit('intpd_mpolel: invalid contraction/shell',0, 175 & BASIS_ERR) 176c 177*** set defNxyz such that it can handle the maximum multi-pole 178c 179 lpole = lval/4 + 1 180 lpole = lpole + 1 ! for derivative 181 call defNxyz(lpole) 182c 183 ucont = (sf_ibs_cn2ucn(ish,ibas)) 184 Li = infbs_cont(CONT_TYPE ,ucont,ibas) 185 i_prim = infbs_cont(CONT_NPRIM,ucont,ibas) 186 i_gen = infbs_cont(CONT_NGEN ,ucont,ibas) 187 i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 188 i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas) 189 i_cent = (sf_ibs_cn2ce(ish,ibas)) 190 i_geom = ibs_geom(ibas) 191c 192 ucont = (sf_ibs_cn2ucn(jsh,jbas)) 193 Lj = infbs_cont(CONT_TYPE ,ucont,jbas) 194 j_prim = infbs_cont(CONT_NPRIM,ucont,jbas) 195 j_gen = infbs_cont(CONT_NGEN ,ucont,jbas) 196 j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas) 197 j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas) 198 j_cent = (sf_ibs_cn2ce(jsh,jbas)) 199 j_geom = ibs_geom(jbas) 200c 201 any_spherical = bas_spherical(ibas).or.bas_spherical(jbas) 202c 203 if (i_geom.ne.j_geom.and.WarnP.eq.0) then 204 write(luout,*)'intpd_mpolel: WARNING: possible geometry', 205 & ' inconsistency' 206 write(luout,*)'i_basis geometry handle:',i_geom 207 write(luout,*)'j_basis geometry handle:',j_geom 208 WarnP = 1 209 endif 210c 211 call intp_txyz(j_cent,j_geom,R,xyz_new) 212c 213 if (i_gen.gt.1 .or. j_gen.gt.1) then 214 write(luout,*) 215 & ' hf3ois does not handle general contractions yet' 216 call errquit('intpd_mpolel: general contraction error ',911, 217 & BASIS_ERR) 218 endif 219c 220 l_int = (lval+1)*(lval+2)/2 221 ij_int = int_nint_cart(i_basis, ish, j_basis, jsh, 0,0, 0,0) 222 num_int = l_int*ij_int 223 num_intd = num_int*3*3 224 if (num_intd.gt.lmpint) then 225 write(luout,*)' intpd_mpolel: lmpint = ',lmpint 226 write(luout,*)' intpd_mpolel: num_intd = ',num_intd 227 call errquit('intpd_mpolel: lmpint too small ',911, 228 & BASIS_ERR) 229 endif 230 if ((lval.eq.0).and.(i_cent.eq.j_cent)) then 231 call dcopy(num_intd,0.0d00,0,MP,1) 232 call ifill(3,0,idatom,1) 233 return 234 endif 235c 236 call hfd3ois( 237 & coords(1,i_cent,i_geom), 238 & dbl_mb(mb_exndcf(i_iexp,ibas)), 239 & dbl_mb(mb_exndcf(i_icfp,ibas)), 240 & i_prim, 1, Li, 241 & xyz_new, 242 & dbl_mb(mb_exndcf(j_iexp,jbas)), 243 & dbl_mb(mb_exndcf(j_icfp,jbas)), 244 & j_prim, 1, Lj, 245 & centerl, 246 & 0.0d00, 247 & 1.0d00, 248 & 1, 1, lval, 249c.....................DryRun 250 & MP,num_int,.false.,scr,lscr) 251 num_mpint = num_intd 252 if (i_cent.eq.j_cent) then 253 idatom(1) = i_cent 254 idatom(2) = 0 255 idatom(3) = -3 256 ioffj = num_int*3 + 1 257 call daxpy(num_int*3,1.0d00,MP(ioffj),1,MP,1) 258 else 259 idatom(1) = i_cent 260 idatom(2) = j_cent 261 idatom(3) = -3 262 endif 263 if (any_spherical) then 264 if (Li.eq.-1) i_gen = 1 265 if (Lj.eq.-1) j_gen = 1 266 call spcart_3cBtran(MP,scr,lscr, 267 & int_nbf_x(Li),int_nbf_s(Li),Li,i_gen,bas_spherical(ibas), 268 & int_nbf_x(Lj),int_nbf_s(Lj),Lj,j_gen,bas_spherical(jbas), 269 & int_nbf_x(lval),int_nbf_x(lval),lval,1,.false., 270 & 9,.false.) 271 endif 272 end 273C> @} 274