1MODULE qs_collocate_density
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 collocate_all_type_1(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
18    TYPE(primitive_type), DIMENSION(:)      :: p
19    REAL                                    :: t1,t2,time
20    INTEGER, DIMENSION(:), ALLOCATABLE      :: r
21
22    INTEGER :: i,j
23
24    ALLOCATE(r(SIZE(p)))
25
26
27    DO i=1,SIZE(rs_grids)
28       rs_grids(i)%r=0.0_dp
29    ENDDO
30
31    CALL CPU_TIME(t1)
32
33    CALL sort_z(rs_grids,p,r)
34    DO i=1,SIZE(p)
35       j=r(i)
36       CALL collocate_pgf_product_rspace(p(j)%la_max,p(j)%zeta,p(j)%la_min, &
37                                         p(j)%lb_max,p(j)%zetb,p(j)%lb_min, &
38                                         p(j)%ra,p(j)%rab,p(j)%rab2, &
39                                         rs_grids(p(j)%igrid),cube_info(p(j)%igrid),l_info,eps_rho_rspace)
40    ENDDO
41
42    CALL CPU_TIME(t2)
43    time=t2-t1
44
45  END SUBROUTINE
46
47  SUBROUTINE sort_z(rs_grids,p,r)
48    TYPE(realspace_grid_type), DIMENSION(:) :: rs_grids
49    TYPE(primitive_type), DIMENSION(:)      :: p
50    INTEGER, DIMENSION(:)      :: r
51    INTEGER, ALLOCATABLE, DIMENSION(:)      :: z
52    INTEGER :: i,n
53    REAL(kind=dp) :: zetp,f,rap(3),rp(3)
54    INTEGER :: cubecenter(3)
55
56    n=SIZE(p)
57    ALLOCATE(z(n))
58    DO i=1,n
59           zetp      = p(i)%zeta + p(i)%zetb
60           f         = p(i)%zetb/zetp
61           rap(:)    = f*p(i)%rab(:)
62           rp(:)     = p(i)%ra(:) + rap(:)
63           cubecenter(:) = MODULO(FLOOR(rp(:)/rs_grids(p(i)%igrid)%dr(:)),rs_grids(p(i)%igrid)%npts(:))
64           z(i)=10000*p(i)%igrid+cubecenter(3)*10+p(i)%la_max+p(i)%lb_max
65    ENDDO
66    CALL sort(z,n,r)
67  END SUBROUTINE
68
69  SUBROUTINE collocate_pgf_product_rspace(la_max,zeta,la_min,&
70                                          lb_max,zetb,lb_min,&
71                                          ra,rab,rab2, &
72                                          rsgrid,cube_info,l_info,&
73                                          eps_rho_rspace)
74
75    INTEGER, INTENT(IN)                      :: la_max
76    REAL(KIND=dp), INTENT(IN)                :: zeta
77    INTEGER, INTENT(IN)                      :: la_min, lb_max
78    REAL(KIND=dp), INTENT(IN)                :: zetb
79    INTEGER, INTENT(IN)                      :: lb_min
80    REAL(KIND=dp), DIMENSION(3), INTENT(IN)  :: ra, rab
81    REAL(KIND=dp), INTENT(IN)                :: rab2
82    TYPE(realspace_grid_type)                :: rsgrid
83    TYPE(cube_info_type), INTENT(IN)         :: cube_info
84    TYPE(l_info_type), INTENT(IN)            :: l_info
85    REAL(KIND=dp), INTENT(IN)                :: eps_rho_rspace
86
87    INTEGER :: cmax, coef_max, gridbounds(2,3), i, ico, icoef, ig, ithread_l, &
88      jco, k, l, la_max_local, la_min_local, lb_max_local, lb_min_local, &
89      length, lx, lx_max, lxa, lxb, lxy, lxy_max, lxyz, lxyz_max, lya, lyb, &
90      lza, lzb, offset, start
91    INTEGER, DIMENSION(3)                    :: cubecenter, lb_cube, ng, &
92                                                ub_cube
93    INTEGER, DIMENSION(:), POINTER           :: ly_max, lz_max, sphere_bounds
94    INTEGER, DIMENSION(:, :), POINTER        :: map
95    LOGICAL                                  :: failure, &
96                                                my_map_consistent
97    REAL(KIND=dp) :: a, b, binomial_k_lxa, binomial_l_lxb, cutoff, f, pg, &
98      prefactor, radius, rpg, ya, yap, yb, ybp, za, zap, zb, zbp, zetp, prefetch
99    REAL(KIND=dp), DIMENSION(3)              :: dr, rap, rb, rbp, roffset, rp
100    REAL(KIND=dp), DIMENSION(:, :, :), &
101      POINTER                                :: grid
102    REAL(KIND=dp), POINTER                   :: alpha_old(:,:), dpy(:,:), &
103                                                dpz(:,:), polx(:,:), &
104                                                poly(:,:), polz(:,:), pzyx(:)
105
106    INTEGER :: lxp,lyp,lzp,lp,iaxis,lxpm,lypm
107    REAL(kind=dp) :: p_ele,ax,ay,az
108    REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: alpha
109    REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xyz
110    REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xyt
111    REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xtt
112
113    REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: pol_z
114    REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: pol_y
115    REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: pol_x
116    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
117
118!   ---------------------------------------------------------------------------
119
120    failure = .FALSE.
121
122    ithread_l=0
123
124    my_map_consistent=.TRUE.
125
126
127    zetp      = zeta + zetb
128    f         = zetb/zetp
129    rap(:)    = f*rab(:)
130    rbp(:)    = rap(:) - rab(:)
131    rp(:)     = ra(:) + rap(:)
132    rb(:)     = ra(:)+rab(:)
133
134    cutoff    = 1.0_dp
135    prefactor = EXP(-zeta*f*rab2)
136    radius=exp_radius_very_extended(la_min,la_max,lb_min,lb_max,ra=ra,rb=rb,rp=rp,&
137                                       zetp=zetp,eps=eps_rho_rspace,&
138                                       prefactor=prefactor,cutoff=cutoff)
139    prefactor = EXP(-zeta*f*rab2)
140
141    IF (radius .EQ. 0.0_dp ) THEN
142      RETURN
143    END IF
144
145        la_max_local=la_max
146        la_min_local=la_min
147        lb_max_local=lb_max
148        lb_min_local=lb_min
149
150    coef_max=la_max_local+lb_max_local+1
151
152!   *** properties of the grid ***
153    dr(:) = rsgrid%dr(:)
154    ng(:) = rsgrid%npts(:)
155
156! WARNING: this resets the lower bounds of grid to 1 (unlike grid => rsgrid%r)
157    grid => rsgrid%r(:,:,:)
158    gridbounds(1,1)=LBOUND(GRID,1)
159    gridbounds(2,1)=UBOUND(GRID,1)
160    gridbounds(1,2)=LBOUND(GRID,2)
161    gridbounds(2,2)=UBOUND(GRID,2)
162    gridbounds(1,3)=LBOUND(GRID,3)
163    gridbounds(2,3)=UBOUND(GRID,3)
164
165!   *** get the sub grid properties for the given radius ***
166    CALL return_cube(cube_info,radius,lb_cube,ub_cube,sphere_bounds)
167
168!   *** get the l_info logic and arrays ***
169    CALL return_l_info(l_info,la_min_local,la_max_local,lb_min_local,lb_max_local,&
170                       ithread_l,lx_max, lxy_max,lxyz_max,ly_max,lz_max, &
171                       map,polx,poly,polz,dpy,dpz,alpha_old,pzyx,cmax)
172
173!   *** position of the gaussian product
174!
175!   this is the actual definition of the position on the grid
176!   i.e. a point rp(:) gets here grid coordinates
177!   MODULO(rp(:)/dr(:),ng(:))+1
178!   hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid ((1,1,1) on grid)
179!
180
181    cubecenter(:) = FLOOR(rp(:)/dr(:))
182    roffset(:)    = rp(:) - REAL(cubecenter(:),dp)*dr(:)
183!   *** a mapping so that the ig corresponds to the right grid point
184    DO i=1,3
185      IF ( rsgrid % perd ( i ) == 1 ) THEN
186        start=lb_cube(i)
187        DO
188         offset=MODULO(cubecenter(i)+start,ng(i))+1-start
189         length=MIN(ub_cube(i),ng(i)-offset)-start
190         DO ig=start,start+length
191            map(ig,i) = ig+offset
192         END DO
193         IF (start+length.GE.ub_cube(i)) EXIT
194         start=start+length+1
195        END DO
196      ELSE
197        ! this takes partial grid + border regions into account
198        offset=MODULO(cubecenter(i),ng(i))+rsgrid%lb(i)
199        offset=offset-rsgrid%lb_local(i)+1
200        DO ig=lb_cube(i),ub_cube(i)
201           map(ig,i) = ig+offset
202        END DO
203      END IF
204    ENDDO
205
206    INCLUDE 'prep.f90'
207
208    ALLOCATE(coef_xyt(((lp+1)*(lp+2))/2))
209    ALLOCATE(coef_xtt(0:lp))
210
211
212!
213!   compute P_{lxp,lyp,lzp} given P_{lxa,lya,lza,lxb,lyb,lzb} and alpha(ls,lxa,lxb,1)
214!   use a three step procedure
215!   we don't store zeros, so counting is done using lxyz,lxy in order to have contiguous memory access in collocate_fast.F
216!
217       lxyz=0
218       DO lzp=0,lp
219       DO lyp=0,lp-lzp
220       DO lxp=0,lp-lzp-lyp
221           lxyz=lxyz+1
222           coef_xyz(lxyz)=0.0_dp
223       ENDDO
224       ENDDO
225       ENDDO
226       DO lzb=0,lb_max_local
227       DO lza=0,la_max_local
228          lxy=0
229          DO lyp=0,lp-lza-lzb
230             DO lxp=0,lp-lza-lzb-lyp
231                lxy=lxy+1
232                coef_xyt(lxy)=0.0_dp
233             ENDDO
234             lxy=lxy+lza+lzb
235          ENDDO
236          DO lyb=0,lb_max_local-lzb
237          DO lya=0,la_max_local-lza
238             lxpm=(lb_max_local-lzb-lyb)+(la_max_local-lza-lya)
239             coef_xtt(0:lxpm)=0.0_dp
240             DO lxb=MAX(lb_min_local-lzb-lyb,0),lb_max_local-lzb-lyb
241             DO lxa=MAX(la_min_local-lza-lya,0),la_max_local-lza-lya
242                ico=coset(lxa,lya,lza)
243                jco=coset(lxb,lyb,lzb)
244                p_ele=prefactor !pab(ico,jco)
245                DO lxp=0,lxa+lxb
246                   coef_xtt(lxp)=coef_xtt(lxp)+p_ele*alpha(lxp,lxa,lxb,1)
247                ENDDO
248             ENDDO
249             ENDDO
250             lxy=0
251             DO lyp=0,lya+lyb
252                DO lxp=0,lp-lza-lzb-lya-lyb
253                  lxy=lxy+1
254                  coef_xyt(lxy)=coef_xyt(lxy)+alpha(lyp,lya,lyb,2)*coef_xtt(lxp)
255                ENDDO
256                lxy=lxy+lza+lzb+lya+lyb-lyp
257             ENDDO
258          ENDDO
259          ENDDO
260          lxyz=0
261          DO lzp=0,lza+lzb
262             lxy=0
263             DO lyp=0,lp-lza-lzb
264                DO lxp=0,lp-lza-lzb-lyp
265                       lxy=lxy+1 ; lxyz=lxyz+1
266                       coef_xyz(lxyz)=coef_xyz(lxyz)+alpha(lzp,lza,lzb,3)*coef_xyt(lxy)
267                ENDDO
268                lxy=lxy+lza+lzb ; lxyz=lxyz+lza+lzb-lzp
269             ENDDO
270             DO lyp=lp-lza-lzb+1,lp-lzp
271                DO lxp=0,lp-lyp-lzp
272                   lxyz=lxyz+1
273                ENDDO
274             ENDDO
275          ENDDO
276       ENDDO
277       ENDDO
278
279    INCLUDE 'call_collocate.f90'
280
281    CALL finish_l_info(polx,poly,polz,pzyx)
282
283  END SUBROUTINE collocate_pgf_product_rspace
284END MODULE
285