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