1 subroutine dim_elfcon(basis,geom,g_dens,points,npt,elfval,nder) 2c 3c $Id: hnd_elfcon.F 19707 2010-10-29 17:59:36Z d3y133 $ 4c 5c This routine calculates the electronic contribution of the 6c electronic integral defined by nder for a given density at 7c the grid points defined in points. 8c 9c It returns an array (max(nder*3,1),npts) which holds all 10c max(nder*3,1) components for each grid point 11c 12 implicit none 13#include "nwc_const.fh" 14#include "errquit.fh" 15#include "global.fh" 16#include "bas.fh" 17#include "mafdecls.fh" 18#include "geom.fh" 19#include "stdio.fh" 20#include "msgids.fh" 21c 22 integer basis ! [input] basis set 23 integer geom ! [input] geometry 24 integer g_dens ! [input] GA with density 25 integer npt ! [input] number of coord points 26 integer nder ! [input] electronic integral type 27 double precision points(3,npt) ! [input] coordinates for points 28 double precision elfval(*) ! [output] efg values for each coord 29c 30 integer ishell, jshell, ijshell, nshell, nbf_max, me, nproc 31 integer ilo, ihi, jlo, jhi, idim, jdim, nint 32 integer l_dens, k_dens, l_scr, k_scr, l_buf, k_buf 33 integer maxbuf, maxscr, i 34c 35 me = ga_nodeid() 36 nproc = ga_nnodes() 37c 38c ----- calculate buffer and scratch space ----- 39c buffer = (lmax*(lmax+1)/2)^2 * (max(nder*3,1) * ngridpoints 40c scratch = see hnd_elfder wrapper routine 41c 42 call int_init_1eelec(maxbuf,maxscr,basis,nder,npt) 43c 44 if (.not. bas_geom(basis, geom)) call errquit 45 $ ('hnd_elfcon: bad basis', 555, BASIS_ERR) 46 if (.not. bas_numcont(basis, nshell)) call errquit 47 $ ('hnd_elfcon: bas_numcont failed for basis', basis, BASIS_ERR) 48 if (.not. bas_nbf_cn_max(basis,nbf_max)) call errquit 49 & ('hnd_elfcon: bas_nbf_cn_max failed',555, BASIS_ERR) 50c 51 if (.not. ma_push_get(mt_dbl,nbf_max*nbf_max,'dens patch',l_dens, 52 & k_dens)) call errquit('hnd_elfcon: ma 1 failed',911,MA_ERR) 53 if (.not. ma_push_get(mt_dbl,maxscr,'scratch',l_scr,k_scr)) 54 & call errquit('hnd_elfcon: ma 2 failed',911,MA_ERR) 55 if (.not. ma_push_get(mt_dbl,maxbuf,'int buf',l_buf,k_buf)) 56 & call errquit('hnd_elfcon: ma 3 failed',911,MA_ERR) 57c 58c Zero elfval result array 59c 60 call dcopy(max(nder*3,1)*npt,0.0d0,0,elfval,1) 61c 62c ----- calculate electronic integral component(s) at all points ----- 63c 64 65 ijshell = 0 66 do ishell = 1, nshell 67c 68c get basis info 69c 70 if (.not. bas_cn2bfr(basis, ishell, ilo, ihi)) call errquit 71 & ('hnd_elfcon: bas_cn2bfr failed for basis',basis,BASIS_ERR) 72 idim = ihi - ilo + 1 73 74 do jshell = 1, nshell 75 ijshell = ijshell + 1 76 if (mod(ijshell,nproc) .eq. me) then 77c 78c get basis info 79c 80 if (.not. bas_cn2bfr(basis, jshell, jlo, jhi)) call 81 & errquit('hnd_elfcon: bas_cn2bfr',basis,BASIS_ERR) 82 jdim = jhi - jlo + 1 83 nint = idim * jdim 84c 85c Get the density patch, make the integrals and contract 86c 87 call ga_get(g_dens, ilo, ihi, jlo, jhi, 88 $ dbl_mb(k_dens), idim) 89c 90 call dim_1eelec(basis,ishell,basis,jshell,maxscr, 91 & dbl_mb(k_scr),nint,dbl_mb(k_buf), 92 & nder,points,npt) 93c 94 call dim_multi_reduce(dbl_mb(k_buf),dbl_mb(k_dens), 95 & elfval,idim,jdim,npt*(max(nder*3,1))) 96 end if ! mod parallel loop 97 end do ! jshell 98 end do ! ishell 99c 100c Collect components from all the nodes for all points 101c 102 call ga_sync() 103 call ga_dgop(msg_efgs_col,elfval,npt*(max(nder*3,1)),'+') 104c 105c Clean up MA data blocks 106c 107 if (.not.ma_pop_stack(l_buf)) call errquit 108 & ('hnd_elfcon, ma_pop_stack of l_buf failed',911,MA_ERR) 109 if (.not.ma_pop_stack(l_scr)) call errquit 110 & ('hnd_elfcon, ma_pop_stack of l_scr failed',911,MA_ERR) 111 if (.not.ma_pop_stack(l_dens)) call errquit 112 & ('hnd_elfcon, ma_pop_stack of l_dens failed',911,MA_ERR) 113 return 114 end 115c 116 subroutine dim_multi_reduce(mblock,block,rblock,idim,jdim,nblock) 117c 118 implicit none 119#include "stdio.fh" 120 integer idim,jdim,nblock 121 double precision mblock(idim,jdim,nblock), block(idim,jdim) 122 double precision rblock(nblock) 123c 124 integer iblock,i,j 125c 126c write(luout,*) "integral patch:" 127c write(luout,*) mblock 128 do iblock = 1, nblock 129 do i = 1, idim 130 do j = 1, jdim 131 rblock(iblock)=rblock(iblock)+mblock(i,j,iblock)* 132 & block(i,j) 133 enddo 134 enddo 135 enddo 136c 137 return 138 end 139