1 subroutine int_giao_1ega(ibas,jbas,g,integ_type,xyzpt,nat, 2 & oskel) 3C$Id$ 4 implicit none 5#include "errquit.fh" 6#include "mafdecls.fh" 7#include "global.fh" 8#include "rtdb.fh" 9#include "inp.fh" 10#include "apiP.fh" 11#include "bas.fh" 12#include "sym.fh" 13#include "geom.fh" 14c 15c Compute the desired type of 1e GIAO integrals 16c and ADD them into the given global array. 17c This version computes the full square of integrals and should work 18c OK even if ibas != jbas. 19c 20c Oskel indicates that the skeleton (petite-list symmetry) matrix should be 21c built ... requires that ibas = jbas. 22c 23c arguments 24c 25 integer ibas, jbas ! [input] bra and ket basis sets 26 integer g ! [output] GA handle to array, one for each field direction (if needed) 27 character*(*) integ_type ! [input] Name of integrals to compute 28 logical oskel ! [input] If true generate symmetry unique list 29 integer nat ! [input] number of atoms for which we get the integrals (if needed) 30 double precision xyzpt(3,*) ! [input] coordinates of requested atoms (if needed) 31c 32c local variables 33c 34 integer mone, two 35 parameter (mone=-1, two=2) 36 integer nshell_i, nshell_j 37 integer ishell, jshell, iproc, nproc, mem1, max1e 38 integer ijshell, ilo, ihi, jlo, jhi, ilen, jlen 39 integer l_buf, l_scr 40 integer k_buf, k_scr 41 integer alo(3), ahi(3), ld(2) 42 integer type 43 logical odoit 44 double precision q2 45 integer nblocks 46c 47 logical odbug 48 logical osome 49c 50 integer irtdb 51 integer int_get_rtdb 52 external int_get_rtdb 53 integer iextbq 54 integer nbq,nextbq,ncosbq,ntmp 55 integer l_qbq,l_cbq 56 integer k_qbq,k_cbq 57 integer k_qextbq,k_cextbq 58c 59 odbug=.false. 60 osome=.false. 61 osome=osome.or.odbug 62 odbug=odbug.and.(ga_nodeid().eq.0) 63 osome=osome.and.(ga_nodeid().eq.0) 64 if(osome) then 65 write(6,*) 'in -int_giao_1ega- ... integ_type = ', 66 $ integ_type,ga_nodeid() 67 call util_flush(6) 68 endif 69c 70 call ga_sync() 71c 72 if (oskel) then 73 if (ibas.ne.jbas) call errquit 74 $ ('int_giao_1ega: use of symmetry requires ibas=jbas', ibas, 75 & BASIS_ERR) 76 end if 77c 78 if (inp_compare(.false., integ_type, 's10')) then 79 type = 1 80 nblocks = 3 81 elseif (inp_compare(.false., integ_type, 'srxRb')) then 82 type = 2 83 nblocks = 3 84 else if (inp_compare(.false., integ_type, 'l10')) then 85 type = 3 86 nblocks = 3 87 else if (inp_compare(.false., integ_type, 'tv10')) then 88 type = 4 89 nblocks = 3 90 else if (inp_compare(.false., integ_type, 'h01')) then 91 type = 5 92 nblocks = 3 * nat 93 else if (inp_compare(.false., integ_type, 'h11 para'))then 94 type = 6 95 nblocks = 9 * nat 96 else if (inp_compare(.false., integ_type, 'h11 dia'))then 97 type = 7 98 nblocks = 9 * nat 99 else if (inp_compare(.false., integ_type, 'h11 all'))then 100 type = 8 101 nblocks = 9 * nat 102 else if (inp_compare(.false., integ_type, 'dso'))then 103 type = 9 104 nblocks = 9 * nat ! nat is number of pairs in this case : nat * (nat-1) 105 else if (inp_compare(.false., integ_type, 'pso'))then 106 type = 10 107 nblocks = 3 * nat 108 else if (inp_compare(.false., integ_type, 'fc'))then 109 type = 11 110 nblocks = nat 111 else if (inp_compare(.false., integ_type, 'sd+fc'))then 112 type = 12 113 nblocks = 6 * nat 114 else if (inp_compare(.false., integ_type, 'velocity'))then 115 type = 13 116 nblocks = 3 117 else if (inp_compare(.false., integ_type, 'angmom'))then 118 type = 14 119 nblocks = 3 120 else if (inp_compare(.false., integ_type, 'bq10')) then 121 type = 15 122 nblocks = 3 123 else 124 write(6,*) ' integ_type = ', integ_type,ga_nodeid() 125 call errquit('int_giao_1ega: unknown integ_type', 0, INT_ERR) 126 end if 127c 128c Get info about the basis sets 129c 130 if (.not. bas_numcont(ibas, nshell_i)) call errquit 131 $ ('rhf_fock_1e: bas_numcont failed for ibas', ibas, 132 & BASIS_ERR) 133 if (.not. bas_numcont(jbas, nshell_j)) call errquit 134 $ ('rhf_fock_1e: bas_numcont failed for jbas', jbas, 135 & BASIS_ERR) 136c 137c allocate necessary local temporary arrays on the stack 138c 139c l_buf ... buffer to hold shell block of matrix 140c l_s ... buffer to hold shell block of matrix 141c l_scr ... workspace for integral routines 142c 143c k_* are the offsets corrsponding to the l_* handles 144c 145c 146 if (type.lt.9) then 147 call int_mem_1e(max1e, mem1) 148 max1e = max1e*nblocks 149 elseif (type.eq.9) then 150 call int_init_dso(max1e,mem1,ibas,nat) 151 elseif (type.eq.10) then 152 call int_init_pso(max1e,mem1,ibas,nat) 153 elseif (type.eq.11.or.type.eq.12) then 154 call int_init_1eelec(max1e,mem1,ibas,2,nat) 155 elseif (type.eq.15) then 156 call int_mem_1e(max1e, mem1) 157 max1e = max1e*nblocks 158 else 159 call int_init_dip(max1e,mem1,ibas) 160 max1e = max1e*nblocks 161 endif 162 mem1 = max(mem1,max1e) 163c 164 if(.not.MA_push_get(MT_DBL,max1e,'int_giao_1ega:buf',l_buf,k_buf)) 165 $ call errquit('int_giao_1ega: ma failed', max1e, MA_ERR) 166 if(.not.MA_push_get(MT_DBL, mem1,'int_giao_1ega:scr',l_scr,k_scr)) 167 $ call errquit('int_giao_1ega: ma failed', mem1, MA_ERR) 168c 169c Get bq charges (external bq and cosmo) 170c 171 if (type.eq.15) then 172 nbq=0 173 nextbq = 0 174 ncosbq = 0 175 if(geom_extbq_on()) then 176 nextbq = geom_extbq_ncenter() 177 k_cextbq = geom_extbq_coord() 178 k_qextbq = geom_extbq_charge() 179 nbq = nextbq ! external bq centers 180 end if 181c 182c Get rtdb handle 183c 184 irtdb = int_get_rtdb() ! get rtdb handle 185 if (rtdb_get(irtdb,'cosmo:nefc',mt_int,1,ncosbq)) 186 & nbq = ncosbq ! cosmo bq centers 187c 188 if (nextbq.gt.0.and.ncosbq.gt.0) 189 & nbq = nextbq + ncosbq ! all bq centers 190c 191c Allocate memory for all the bq charges 192c 193 if(.not.MA_push_get(MT_DBL,nbq,'qbq',l_qbq,k_qbq)) 194 & call errquit('int_giao_1ega: ma failed', l_qbq, MA_ERR) 195 call dfill(nbq,0.d0,dbl_mb(k_qbq),1) 196 if(.not.MA_push_get(MT_DBL,3*nbq,'cbq',l_cbq,k_cbq)) 197 $ call errquit('int_giao_1ega: ma failed', l_cbq, MA_ERR) 198 call dfill(3*nbq,0.d0,dbl_mb(k_cbq),1) 199c 200c Assign external bq charges and coordinates 201c 202 if (nextbq.gt.0) then 203 call dcopy(nextbq,dbl_mb(k_qextbq),1,dbl_mb(k_qbq),1) 204 call dcopy(3*nextbq,dbl_mb(k_cextbq),1,dbl_mb(k_cbq),1) 205 end if ! nextbq 206c 207c Get cosmo charges and coordinates 208c 209 if (ncosbq.gt.0) then 210 call bq_fromrtdb(irtdb,'cosmo:efcz','cosmo:efcc', 211 & ntmp,dbl_mb(k_qbq+nextbq),dbl_mb(k_cbq+3*nextbq)) 212 end if ! ncosbq gt 0 213 end if ! type .eq. 15 214c 215c Loop thru shells with static parallel work decomposition 216c 217 iproc = ga_nodeid() 218 nproc = ga_nnodes() 219 ijshell = 0 220 q2 = 1.0d0 221 do jshell = 1, nshell_j 222 do ishell = 1, nshell_i 223c 224 if (mod(ijshell, nproc) .eq. iproc) then 225 odoit = .true. 226 if (oskel) 227 $ odoit = sym_shell_pair(ibas, ishell, jshell, q2) 228c 229 if (odoit) then 230 if (.not. bas_cn2bfr(ibas, ishell, ilo, ihi)) 231 $ call errquit('int_1e_ga: bas_cn2bfr ?', ibas, 232 & BASIS_ERR) 233 if (.not. bas_cn2bfr(jbas, jshell, jlo, jhi)) 234 $ call errquit('int_1e_ga: bas_cn2bfr ?', jbas, 235 & BASIS_ERR) 236c 237 ilen = ihi-ilo+1 238 jlen = jhi-jlo+1 239c 240c Generate the integrals 241c 242 if (type .eq. 1) then ! 3 243 call int_giaos10 (jbas, jshell, ibas, ishell, 244 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf)) 245 else if (type .eq. 2) then ! 3 246 call int_giaos100(jbas, jshell, ibas, ishell, 247 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf)) 248 else if (type .eq. 3) then ! 3 249 call int_giaol10 (jbas, jshell, ibas, ishell, 250 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf)) 251 else if (type .eq. 4) then ! 3 252 call int_giaotv10 (jbas, jshell, ibas, ishell, 253 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf)) 254 else if (type .eq. 5) then ! 3*nat 255 call int_giaoh01 (jbas, jshell, ibas, ishell, 256 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 257 $ xyzpt, nat) 258 else if (type .eq. 6) then ! 9*nat 259 call int_giaoh11 (jbas, jshell, ibas, ishell, 260 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 261 $ xyzpt, nat, .true.,.false.) 262 else if (type .eq. 7) then ! 9*nat 263 call int_giaoh11 (jbas, jshell, ibas, ishell, 264 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 265 $ xyzpt, nat, .false.,.true.) 266 else if (type .eq. 8) then ! 9*nat 267 call int_giaoh11 (jbas, jshell, ibas, ishell, 268 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 269 $ xyzpt, nat, .true.,.true.) 270 else if (type .eq. 9) then ! 9*nat*(nat-1) 271 call int_dso (jbas, jshell, ibas, ishell, 272 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 273 $ xyzpt, nat) 274 else if (type .eq. 10) then ! 3*nat 275 call int_pso (jbas, jshell, ibas, ishell, 276 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 277 $ xyzpt, nat) 278 else if (type .eq. 11) then ! nat 279 call int_1eelec(jbas, jshell, ibas,ishell, 280 & mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 281 & mone,xyzpt,nat) 282 else if (type .eq. 12) then ! 6*nat 283 call int_1eelec(jbas, jshell, ibas,ishell, 284 & mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 285 & two,xyzpt,nat) 286 else if (type .eq. 13) then ! 3 287 call int_veloc(jbas, jshell, ibas, ishell, 288 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf)) 289 else if (type .eq. 14) then ! 3 290 call int_angmom(jbas, jshell, ibas, ishell, 291 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf)) 292 else if (type .eq. 15) then ! 3 293 call int_giaobq10(jbas, jshell, ibas, ishell, 294 $ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf), 295 & dbl_mb(k_qbq),dbl_mb(k_cbq),nbq) 296 else 297 call errquit('int_giao_1ega: invalid type?', type, 298 & GA_ERR) 299 end if 300c 301c Add the integrals into the global array 302c 303 alo(1) = jlo 304 ahi(1) = jhi 305 alo(2) = ilo 306 ahi(2) = ihi 307 alo(3) = 1 308 ahi(3) = nblocks 309 ld(1) = jlen 310 ld(2) = ilen 311 call nga_acc(g,alo,ahi,dbl_mb(k_buf),ld,1.0d0) 312 end if 313 endif 314 ijshell = ijshell + 1 315 end do 316 end do 317c 318c chop stack at first item allocated 319c 320 if (type.eq.15) then ! for bq's 321 if (.not. MA_pop_stack(l_cbq)) call errquit 322 $ ('int_giao_1ega: pop failed', l_cbq, MA_ERR) 323 if (.not. MA_pop_stack(l_qbq)) call errquit 324 $ ('int_giao_1ega: pop failed', l_qbq, MA_ERR) 325 end if ! type 15 326 if (.not. MA_pop_stack(l_scr)) call errquit 327 $ ('int_giao_1ega: pop failed', l_scr, MA_ERR) 328 if (.not. MA_pop_stack(l_buf)) call errquit 329 $ ('int_giao_1ega: pop failed', l_buf, MA_ERR) 330 331 call ga_sync() ! So that no nasty races can result 332c 333 end 334