1MODULE qs_integrate_potential 2 use kinds 3 use util 4 use l_utils 5 use cube_utils 6 use qs_interactions 7 use orbital_pointers 8 use basic_types 9 IMPLICIT NONE 10 11CONTAINS 12 13 SUBROUTINE integrate_all(p,rs_grids,cube_info,l_info,eps_rho_rspace,time) 14 TYPE(realspace_grid_type), DIMENSION(:) :: rs_grids 15 TYPE(cube_info_type), DIMENSION(:) :: cube_info 16 TYPE(l_info_type) :: l_info 17 REAL(KIND=dp) :: eps_rho_rspace,intsum 18 TYPE(primitive_type), DIMENSION(:) :: p 19 REAL :: t1,t2,time 20 INTEGER, DIMENSION(:), ALLOCATABLE :: r 21 22 INTEGER :: i,j,k,l,m 23 24 ALLOCATE(r(SIZE(p))) 25 26 DO i=1,SIZE(rs_grids) 27 DO k=LBOUND(rs_grids(i)%r,3),UBOUND(rs_grids(i)%r,3) 28 DO l=LBOUND(rs_grids(i)%r,2),UBOUND(rs_grids(i)%r,2) 29 DO m=LBOUND(rs_grids(i)%r,1),UBOUND(rs_grids(i)%r,1) 30 rs_grids(i)%r(m,l,k)=1.0_dp/sqrt(1.0_dp+(REAL(m+3*k+5*l,dp))**2) 31 ENDDO 32 ENDDO 33 ENDDO 34 ENDDO 35 36 CALL CPU_TIME(t1) 37 38 DO i=1,SIZE(p) 39 j=i 40 CALL integrate_pgf_product_rspace(p(j)%la_max,p(j)%zeta,p(j)%la_min, & 41 p(j)%lb_max,p(j)%zetb,p(j)%lb_min, & 42 p(j)%ra,p(j)%rab,p(j)%rab2, & 43 rs_grids(p(j)%igrid),cube_info(p(j)%igrid),l_info,eps_rho_rspace,& 44 intsum) 45 p(j)%intsum=intsum 46 ENDDO 47 48 CALL CPU_TIME(t2) 49 time=t2-t1 50 51 END SUBROUTINE 52 53 SUBROUTINE integrate_pgf_product_rspace(la_max,zeta,la_min,& 54 lb_max,zetb,lb_min,& 55 ra,rab,rab2, & 56 rsgrid,cube_info,l_info,& 57 eps_rho_rspace,& 58 intsum) 59 60 INTEGER, INTENT(IN) :: la_max 61 REAL(KIND=dp), INTENT(IN) :: zeta 62 INTEGER, INTENT(IN) :: la_min, lb_max 63 REAL(KIND=dp), INTENT(IN) :: zetb 64 INTEGER, INTENT(IN) :: lb_min 65 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab 66 REAL(KIND=dp), INTENT(IN) :: rab2 67 TYPE(realspace_grid_type) :: rsgrid 68 TYPE(cube_info_type), INTENT(IN) :: cube_info 69 TYPE(l_info_type), INTENT(IN) :: l_info 70 REAL(KIND=dp), INTENT(IN) :: eps_rho_rspace 71 REAL(KIND=dp), INTENT(OUT) :: intsum 72 73 INTEGER :: cmax, coef_max, gridbounds(2,3), i, ico, icoef, ig, ithread_l, & 74 jco, k, l, la_max_local, la_min_local, lb_max_local, lb_min_local, & 75 length, lx, lx_max, lxa, lxb, lxy, lxy_max, lxyz, lxyz_max, lya, lyb, & 76 lza, lzb, offset, start 77 INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, & 78 ub_cube 79 INTEGER, DIMENSION(:), POINTER :: ly_max, lz_max, sphere_bounds 80 INTEGER, DIMENSION(:, :), POINTER :: map 81 LOGICAL :: failure, & 82 my_map_consistent 83 REAL(KIND=dp) :: a, b, binomial_k_lxa, binomial_l_lxb, cutoff, pg, & 84 prefactor, radius, rpg, ya, yap, yb, ybp, za, zap, zb, zbp, zetp, prefetch, & 85 exp_x0, exp_x1, exp_x2, f, ftza, ftzb 86 REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp 87 REAL(KIND=dp), DIMENSION(:, :, :), & 88 POINTER :: grid 89 REAL(KIND=dp), POINTER :: alpha_old(:,:), dpy(:,:), & 90 dpz(:,:), polx(:,:), & 91 poly(:,:), polz(:,:), pzyx(:) 92 93 INTEGER :: lxp,lyp,lzp,lp,iaxis 94 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: alpha 95 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xyz 96 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xyt 97 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xtt 98 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: coef_ttz 99 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: coef_tyz 100 101 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: pol_z 102 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: pol_y 103 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: pol_x 104 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: vab 105 REAL(KIND=dp) :: t_exp_1,t_exp_2,t_exp_min_1,t_exp_min_2,t_exp_plus_1,t_exp_plus_2 106 107 108! --------------------------------------------------------------------------- 109 110 failure = .FALSE. 111 112 ithread_l=0 113 114 my_map_consistent=.TRUE. 115 116 117 zetp = zeta + zetb 118 f = zetb/zetp 119 rap(:) = f*rab(:) 120 rbp(:) = rap(:) - rab(:) 121 rp(:) = ra(:) + rap(:) 122 rb(:) = ra(:)+rab(:) 123 124 cutoff = 1.0_dp 125 prefactor = EXP(-zeta*f*rab2) 126 radius=exp_radius_very_extended(la_min,la_max,lb_min,lb_max,ra=ra,rb=rb,rp=rp,& 127 zetp=zetp,eps=eps_rho_rspace,& 128 prefactor=prefactor,cutoff=cutoff) 129 prefactor = EXP(-zeta*f*rab2) 130 131 IF (radius .EQ. 0.0_dp ) THEN 132 RETURN 133 END IF 134 135 la_max_local=la_max 136 la_min_local=la_min 137 lb_max_local=lb_max 138 lb_min_local=lb_min 139 140 coef_max=la_max_local+lb_max_local+1 141 142! *** properties of the grid *** 143 dr(:) = rsgrid%dr(:) 144 ng(:) = rsgrid%npts(:) 145 146! WARNING: this resets the lower bounds of grid to 1 (unlike grid => rsgrid%r) 147 grid => rsgrid%r(:,:,:) 148 gridbounds(1,1)=LBOUND(GRID,1) 149 gridbounds(2,1)=UBOUND(GRID,1) 150 gridbounds(1,2)=LBOUND(GRID,2) 151 gridbounds(2,2)=UBOUND(GRID,2) 152 gridbounds(1,3)=LBOUND(GRID,3) 153 gridbounds(2,3)=UBOUND(GRID,3) 154 155 CALL return_cube(cube_info,radius,lb_cube,ub_cube,sphere_bounds) 156 CALL return_l_info(l_info,la_min_local,la_max_local,lb_min_local,lb_max_local, & 157 ithread_l,lx_max,lxy_max,lxyz_max,ly_max,lz_max,& 158 map,polx,poly,polz,dpy,dpz,alpha_old,pzyx,cmax) 159 160 dr(:) = rsgrid%dr(:) 161 ng(:) = rsgrid%npts(:) 162 163 grid => rsgrid%r(:,:,:) 164 gridbounds(1,1)=LBOUND(GRID,1) 165 gridbounds(2,1)=UBOUND(GRID,1) 166 gridbounds(1,2)=LBOUND(GRID,2) 167 gridbounds(2,2)=UBOUND(GRID,2) 168 gridbounds(1,3)=LBOUND(GRID,3) 169 gridbounds(2,3)=UBOUND(GRID,3) 170 171 cubecenter(:) = FLOOR(rp(:)/dr(:)) 172 roffset(:) = rp(:) - REAL(cubecenter(:),dp)*dr(:) 173 174 DO i=1,3 175 IF ( rsgrid % perd ( i ) == 1 ) THEN 176 start=lb_cube(i) 177 DO 178 offset=MODULO(cubecenter(i)+start,ng(i))+1-start 179 length=MIN(ub_cube(i),ng(i)-offset)-start 180 DO ig=start,start+length 181 map(ig,i) = ig+offset 182 END DO 183 IF (start+length.GE.ub_cube(i)) EXIT 184 start=start+length+1 185 END DO 186 ELSE 187 ! this takes partial grid + border regions into account 188 offset=MODULO(cubecenter(i),ng(i))+rsgrid%lb(i) 189 offset=offset-rsgrid%lb_local(i)+1 190 DO ig=lb_cube(i),ub_cube(i) 191 map(ig,i) = ig+offset 192 END DO 193 END IF 194 ENDDO 195 196 INCLUDE 'prep.f90' 197 198 INCLUDE 'call_integrate.f90' 199 200! CALL integrate_core_default(grid(1,1,1),coef_xyz(1),pol_x(0,-cmax),pol_y(1,0,-cmax),pol_z(1,0,-cmax), & 201! map(-cmax,1),sphere_bounds(1),lp,cmax,gridbounds(1,1)) 202 203 ! 204 ! compute v_{lxa,lya,lza,lxb,lyb,lzb} given v_{lxp,lyp,lzp} and alpha(ls,lxa,lxb,1) 205 ! use a three step procedure 206 ! 207 208 ALLOCATE(vab(ncoset(la_max_local),ncoset(lb_max_local))) 209 ALLOCATE(coef_ttz(0:la_max_local,0:lb_max_local)) 210 ALLOCATE(coef_tyz(0:la_max_local,0:lb_max_local,0:la_max_local,0:lb_max_local)) 211 212 213 IF (.TRUE.) THEN 214 215 vab=0.0_dp 216 lxyz=0 217 DO lzp=0,lp 218 coef_tyz=0.0_dp 219 DO lyp=0,lp-lzp 220 coef_ttz=0.0_dp 221 DO lxp=0,lp-lzp-lyp 222 lxyz=lxyz+1 223 DO lxb=0,lb_max_local 224 DO lxa=0,la_max_local 225 coef_ttz(lxa,lxb)=coef_ttz(lxa,lxb)+coef_xyz(lxyz)*alpha(lxp,lxa,lxb,1) 226 ENDDO 227 ENDDO 228 229 ENDDO 230 231 DO lyb=0,lb_max_local 232 DO lya=0,la_max_local 233 DO lxb=0,lb_max_local-lyb 234 DO lxa=0,la_max_local-lya 235 coef_tyz(lxa,lxb,lya,lyb)=coef_tyz(lxa,lxb,lya,lyb)+coef_ttz(lxa,lxb)*alpha(lyp,lya,lyb,2) 236 ENDDO 237 ENDDO 238 ENDDO 239 ENDDO 240 241 ENDDO 242 243 DO lzb=0,lb_max_local 244 DO lza=0,la_max_local 245 DO lyb=0,lb_max_local-lzb 246 DO lya=0,la_max_local-lza 247 DO lxb=MAX(lb_min_local-lzb-lyb,0),lb_max_local-lzb-lyb 248 jco=coset(lxb,lyb,lzb) 249 DO lxa=MAX(la_min_local-lza-lya,0),la_max_local-lza-lya 250 ico=coset(lxa,lya,lza) 251 vab(ico,jco)=vab(ico,jco)+coef_tyz(lxa,lxb,lya,lyb)*alpha(lzp,lza,lzb,3) 252 ENDDO 253 ENDDO 254 ENDDO 255 ENDDO 256 ENDDO 257 ENDDO 258 259 ENDDO 260 261 ELSE 262 263 DO lzb=0,lb_max_local 264 DO lza=0,la_max_local 265 DO lyb=0,lb_max_local-lzb 266 DO lya=0,la_max_local-lza 267 DO lxb=MAX(lb_min_local-lzb-lyb,0),lb_max_local-lzb-lyb 268 DO lxa=MAX(la_min_local-lza-lya,0),la_max_local-lza-lya 269 ico=coset(lxa,lya,lza) 270 jco=coset(lxb,lyb,lzb) 271 vab(ico,jco)=0.0_dp 272 lxyz=0 273 DO lzp=0,lp 274 DO lyp=0,lp-lzp 275 DO lxp=0,lp-lzp-lyp 276 lxyz=lxyz+1 277 vab(ico,jco)=vab(ico,jco)+coef_xyz(lxyz)*alpha(lzp,lza,lzb,3)*alpha(lyp,lya,lyb,2)*alpha(lxp,lxa,lxb,1) 278 ENDDO 279 ENDDO 280 ENDDO 281 ENDDO 282 ENDDO 283 ENDDO 284 ENDDO 285 ENDDO 286 ENDDO 287 288 ENDIF 289 290 ! compute a kind of forces checksum 291 intsum=0.0_dp 292 lxyz=0 293 DO lxa=0,la_max 294 DO lxb=0,lb_max 295 DO lya=0,la_max-lxa 296 DO lyb=0,lb_max-lxb 297 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 298 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 299 ico=coset(lxa,lya,lza) 300 jco=coset(lxb,lyb,lzb) 301 intsum=intsum+vab(ico,jco)*(coset(lxa,lya,lza)+1000*coset(lxb,lyb,lzb)) 302 ENDDO 303 ENDDO 304 ENDDO 305 ENDDO 306 ENDDO 307 ENDDO 308 309 CALL finish_l_info(polx,poly,polz,pzyx) 310 311 END SUBROUTINE integrate_pgf_product_rspace 312END MODULE qs_integrate_potential 313