Skip to content

Commit

Permalink
grid: Move collocate_gauge_ortho into qs_linres_current.F
Browse files Browse the repository at this point in the history
  • Loading branch information
@oschuett
oschuett committed Feb 14, 2020
1 parent e1bf4e7 commit e4387f1
Show file tree
Hide file tree
Showing 2 changed files with 163 additions and 127 deletions.
123 changes: 3 additions & 120 deletions src/grid/grid_collocate.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,6 @@ MODULE grid_collocate
!> \param radius ...
!> \param idir ...
!> \param ir ...
!> \param rsgauge ...
!> \param rsbuf ...
!> \param use_subpatch ...
!> \param subpatch_pattern ...
!> \param lmax_global Maximum possible value of lmax used to dimension arrays
Expand All @@ -82,7 +80,7 @@ SUBROUTINE collocate_pgf_product(la_max, zeta, la_min, &
ra, rab, scale, pab, o1, o2, &
rsgrid, cell, cube_info, &
ga_gb_function, &
radius, idir, ir, rsgauge, rsbuf, &
radius, idir, ir, &
use_subpatch, subpatch_pattern, &
lmax_global)

Expand All @@ -101,7 +99,6 @@ SUBROUTINE collocate_pgf_product(la_max, zeta, la_min, &
REAL(KIND=dp), INTENT(IN) :: radius
INTEGER, INTENT(IN) :: ga_gb_function
INTEGER, INTENT(IN), OPTIONAL :: idir, ir
TYPE(realspace_grid_type), POINTER, OPTIONAL :: rsgauge, rsbuf
LOGICAL, OPTIONAL :: use_subpatch
INTEGER(KIND=int_8), OPTIONAL, INTENT(IN):: subpatch_pattern
INTEGER, INTENT(IN) :: lmax_global
Expand Down Expand Up @@ -136,8 +133,6 @@ SUBROUTINE collocate_pgf_product(la_max, zeta, la_min, &
REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pol_x
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

REAL(kind=dp), POINTER, DIMENSION(:, :, :) :: grid_tmp, gauge

subpatch_collocate = .FALSE.

IF (PRESENT(use_subpatch)) THEN
Expand Down Expand Up @@ -406,12 +401,6 @@ SUBROUTINE collocate_pgf_product(la_max, zeta, la_min, &
gridbounds(1, 3) = LBOUND(GRID, 3)
gridbounds(2, 3) = UBOUND(GRID, 3)
IF (PRESENT(ir) .AND. PRESENT(rsgauge) .AND. PRESENT(rsbuf)) THEN
grid => rsbuf%r(:, :, :)
grid_tmp => rsgrid%r(:, :, :)
gauge => rsgauge%r(:, :, :)
ENDIF
! *** initialise the coefficient matrix, we transform the sum
!
! sum_{lxa,lya,lza,lxb,lyb,lzb} P_{lxa,lya,lza,lxb,lyb,lzb} *
Expand Down Expand Up @@ -595,13 +584,9 @@ SUBROUTINE collocate_ortho()
ALLOCATE (pol_y(1:2, 0:lp, -cmax:0))
ALLOCATE (pol_x(0:lp, -cmax:cmax))
IF (PRESENT(ir) .AND. PRESENT(rsgauge)) CALL collocate_ortho_set_to_0()
#include "prep.f90"
#include "call_collocate.f90"
IF (PRESENT(ir) .AND. PRESENT(rsgauge)) CALL collocate_gauge_ortho()
! deallocation needed to pass around a pgi bug..
DEALLOCATE (pol_z)
DEALLOCATE (pol_y)
Expand All @@ -611,106 +596,7 @@ SUBROUTINE collocate_ortho()
END SUBROUTINE collocate_ortho
! **************************************************************************************************
!> \brief ...
! **************************************************************************************************
SUBROUTINE collocate_gauge_ortho()
INTEGER :: i, igmax, igmin, j, j2, jg, jg2, jgmin, &
k, k2, kg, kg2, kgmin, sci
REAL(KIND=dp) :: point(3, 4), res(4), x, y, y2, z, z2
! notice we're in the ortho case

dr(1) = rsgrid%desc%dh(1, 1)
dr(2) = rsgrid%desc%dh(2, 2)
dr(3) = rsgrid%desc%dh(3, 3)
!
sci = 1
kgmin = sphere_bounds(sci)
sci = sci + 1
DO kg = kgmin, 0
kg2 = 1 - kg
k = map(kg, 3)
k2 = map(kg2, 3)
jgmin = sphere_bounds(sci)
sci = sci + 1
z = (REAL(kg, dp) + REAL(cubecenter(3), dp))*dr(3)
z2 = (REAL(kg2, dp) + REAL(cubecenter(3), dp))*dr(3)
DO jg = jgmin, 0
jg2 = 1 - jg
j = map(jg, 2)
j2 = map(jg2, 2)
igmin = sphere_bounds(sci)
sci = sci + 1
igmax = 1 - igmin
y = (REAL(jg, dp) + REAL(cubecenter(2), dp))*dr(2)
y2 = (REAL(jg2, dp) + REAL(cubecenter(2), dp))*dr(2)
DO ig = igmin, igmax
i = map(ig, 1)
x = (REAL(ig, dp) + REAL(cubecenter(1), dp))*dr(1)
point(1, 1) = x; point(2, 1) = y; point(3, 1) = z
point(1, 2) = x; point(2, 2) = y2; point(3, 2) = z
point(1, 3) = x; point(2, 3) = y; point(3, 3) = z2
point(1, 4) = x; point(2, 4) = y2; point(3, 4) = z2
!
res(1) = (point(ir, 1) - rb(ir)) - gauge(i, j, k)
res(2) = (point(ir, 2) - rb(ir)) - gauge(i, j2, k)
res(3) = (point(ir, 3) - rb(ir)) - gauge(i, j, k2)
res(4) = (point(ir, 4) - rb(ir)) - gauge(i, j2, k2)
!
grid_tmp(i, j, k) = grid_tmp(i, j, k) + grid(i, j, k)*res(1)
grid_tmp(i, j2, k) = grid_tmp(i, j2, k) + grid(i, j2, k)*res(2)
grid_tmp(i, j, k2) = grid_tmp(i, j, k2) + grid(i, j, k2)*res(3)
grid_tmp(i, j2, k2) = grid_tmp(i, j2, k2) + grid(i, j2, k2)*res(4)
ENDDO
ENDDO
ENDDO
END SUBROUTINE collocate_gauge_ortho

! **************************************************************************************************
!> \brief ...
! **************************************************************************************************
SUBROUTINE collocate_ortho_set_to_0()
INTEGER :: i, igmax, igmin, j, j2, jg, jg2, jgmin, &
k, k2, kg, kg2, kgmin, sci

!

dr(1) = rsgrid%desc%dh(1, 1)
dr(2) = rsgrid%desc%dh(2, 2)
dr(3) = rsgrid%desc%dh(3, 3)
!
sci = 1
kgmin = sphere_bounds(sci)
sci = sci + 1
DO kg = kgmin, 0
kg2 = 1 - kg
k = map(kg, 3)
k2 = map(kg2, 3)
jgmin = sphere_bounds(sci)
sci = sci + 1
DO jg = jgmin, 0
jg2 = 1 - jg
j = map(jg, 2)
j2 = map(jg2, 2)
igmin = sphere_bounds(sci)
sci = sci + 1
igmax = 1 - igmin
DO ig = igmin, igmax
i = map(ig, 1)
grid(i, j, k) = 0.0_dp
grid(i, j2, k) = 0.0_dp
grid(i, j, k2) = 0.0_dp
grid(i, j2, k2) = 0.0_dp
ENDDO
ENDDO
ENDDO
END SUBROUTINE collocate_ortho_set_to_0

!
! this is a general 'optimized' routine to do the collocation
!
! **************************************************************************************************
!> \brief ...
!> \brief this is a general 'optimized' routine to do the collocation
! **************************************************************************************************
SUBROUTINE collocate_general_opt()
Expand Down Expand Up @@ -977,11 +863,8 @@ SUBROUTINE collocate_general_wings()
END SUBROUTINE
!
! this is a general 'reference' routine to do the collocation
!
! **************************************************************************************************
!> \brief ...
!> \brief this is a general 'reference' routine to do the collocation
! **************************************************************************************************
SUBROUTINE collocate_general()
Expand Down
167 changes: 160 additions & 7 deletions src/qs_linres_current.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,9 @@ MODULE qs_linres_current
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
USE cube_utils, ONLY: cube_info_type
USE cube_utils, ONLY: compute_cube_center,&
cube_info_type,&
return_cube
USE dbcsr_api, ONLY: &
dbcsr_add_block_node, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, &
Expand Down Expand Up @@ -1149,30 +1151,43 @@ SUBROUTINE calculate_jrho_resp(mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, iB, idir
! f_{ab} = g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} -
! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b}
scale = 1.0_dp
current_env%rs_buf(igrid_level)%rs_grid%r(:, :, :) = 0.0_dp
CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
rs_current(igrid_level)%rs_grid, cell, cube_info(igrid_level), &
cell=cell, cube_info=cube_info(igrid_level), &
radius=radius, &
ga_gb_function=GRID_FUNC_ADBmDAB, &
idir=idir, ir=iiiB, &
rsgauge=rs_gauge(iiiB)%rs(igrid_level)%rs_grid, &
rsbuf=current_env%rs_buf(igrid_level)%rs_grid, &
rsgrid=current_env%rs_buf(igrid_level)%rs_grid, &
lmax_global=lmax_global)
CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level)%rs_grid, &
rsbuf=rs_current(igrid_level)%rs_grid, &
rsgauge=rs_gauge(iiiB)%rs(igrid_level)%rs_grid, &
cube_info=cube_info(igrid_level), radius=radius, &
zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
ra=ra, rab=rab, ir=iiiB)

! here the decontracted mat_jp_riii{ab} is multiplied by
! f_{ab} = -g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} +
! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b}
scale = -1.0_dp
current_env%rs_buf(igrid_level)%rs_grid%r(:, :, :) = 0.0_dp
CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
rs_current(igrid_level)%rs_grid, cell, cube_info(igrid_level), &
cell=cell, cube_info=cube_info(igrid_level), &
radius=radius, &
ga_gb_function=GRID_FUNC_ADBmDAB, &
idir=idir, ir=iiB, &
rsgauge=rs_gauge(iiB)%rs(igrid_level)%rs_grid, &
rsbuf=current_env%rs_buf(igrid_level)%rs_grid, &
rsgrid=current_env%rs_buf(igrid_level)%rs_grid, &
lmax_global=lmax_global)
CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level)%rs_grid, &
rsbuf=rs_current(igrid_level)%rs_grid, &
rsgauge=rs_gauge(iiB)%rs(igrid_level)%rs_grid, &
cube_info=cube_info(igrid_level), radius=radius, &
zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
ra=ra, rab=rab, ir=iiB)
ELSE
! here the decontracted mat_jp_rii{ab} is multiplied by
! f_{ab} = g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} -
Expand Down Expand Up @@ -1240,6 +1255,144 @@ SUBROUTINE calculate_jrho_resp(mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, iB, idir

END SUBROUTINE calculate_jrho_resp

! **************************************************************************************************
!> \brief ...
!> \param rsgrid ...
!> \param rsbuf ...
!> \param rsgauge ...
!> \param cube_info ...
!> \param radius ...
!> \param ra ...
!> \param rab ...
!> \param zeta ...
!> \param zetb ...
!> \param ir ...
! **************************************************************************************************
SUBROUTINE collocate_gauge_ortho(rsgrid, rsbuf, rsgauge, cube_info, radius, ra, rab, zeta, zetb, ir)
TYPE(realspace_grid_type) :: rsgrid, rsbuf, rsgauge
TYPE(cube_info_type), INTENT(IN) :: cube_info
REAL(KIND=dp), INTENT(IN) :: radius
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab
REAL(KIND=dp), INTENT(IN) :: zeta, zetb
INTEGER, INTENT(IN) :: ir

INTEGER :: cmax, i, ig, igmax, igmin, j, j2, jg, &
jg2, jgmin, k, k2, kg, kg2, kgmin, &
length, offset, sci, start
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map
INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, ub_cube
INTEGER, DIMENSION(:), POINTER :: sphere_bounds
REAL(KIND=dp) :: f, point(3, 4), res(4), x, y, y2, z, z2, &
zetp
REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp
REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gauge, grid, grid_buf

CPASSERT(rsgrid%desc%orthorhombic)
NULLIFY (sphere_bounds)

grid => rsgrid%r(:, :, :)
grid_buf => rsbuf%r(:, :, :)
gauge => rsgauge%r(:, :, :)

! *** center of gaussians and their product
zetp = zeta + zetb
f = zetb/zetp
rap(:) = f*rab(:)
rbp(:) = rap(:) - rab(:)
rp(:) = ra(:) + rap(:)
rb(:) = ra(:) + rab(:)

! *** properties of the grid ***
ng(:) = rsgrid%desc%npts(:)
dr(1) = rsgrid%desc%dh(1, 1)
dr(2) = rsgrid%desc%dh(2, 2)
dr(3) = rsgrid%desc%dh(3, 3)

! *** get the sub grid properties for the given radius ***
CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
cmax = MAXVAL(ub_cube)

! *** position of the gaussian product
!
! this is the actual definition of the position on the grid
! i.e. a point rp(:) gets here grid coordinates
! MODULO(rp(:)/dr(:),ng(:))+1
! hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid ((1,1,1) on grid)
!
ALLOCATE (map(-cmax:cmax, 3))
map(:, :) = -1
CALL compute_cube_center(cubecenter, rsgrid%desc, zeta, zetb, ra, rab)
roffset(:) = rp(:) - REAL(cubecenter(:), dp)*dr(:)

! *** a mapping so that the ig corresponds to the right grid point
DO i = 1, 3
IF (rsgrid%desc%perd(i) == 1) THEN
start = lb_cube(i)
DO
offset = MODULO(cubecenter(i) + start, ng(i)) + 1 - start
length = MIN(ub_cube(i), ng(i) - offset) - start
DO ig = start, start + length
map(ig, i) = ig + offset
END DO
IF (start + length .GE. ub_cube(i)) EXIT
start = start + length + 1
END DO
ELSE
! this takes partial grid + border regions into account
offset = MODULO(cubecenter(i) + lb_cube(i) + rsgrid%desc%lb(i) - rsgrid%lb_local(i), ng(i)) + 1 - lb_cube(i)
! check for out of bounds
IF (ub_cube(i) + offset > UBOUND(grid, i) .OR. lb_cube(i) + offset < LBOUND(grid, i)) THEN
CPASSERT(.FALSE.)
ENDIF
DO ig = lb_cube(i), ub_cube(i)
map(ig, i) = ig + offset
END DO
END IF
ENDDO

! *** actually loop over the grid
sci = 1
kgmin = sphere_bounds(sci)
sci = sci + 1
DO kg = kgmin, 0
kg2 = 1 - kg
k = map(kg, 3)
k2 = map(kg2, 3)
jgmin = sphere_bounds(sci)
sci = sci + 1
z = (REAL(kg, dp) + REAL(cubecenter(3), dp))*dr(3)
z2 = (REAL(kg2, dp) + REAL(cubecenter(3), dp))*dr(3)
DO jg = jgmin, 0
jg2 = 1 - jg
j = map(jg, 2)
j2 = map(jg2, 2)
igmin = sphere_bounds(sci)
sci = sci + 1
igmax = 1 - igmin
y = (REAL(jg, dp) + REAL(cubecenter(2), dp))*dr(2)
y2 = (REAL(jg2, dp) + REAL(cubecenter(2), dp))*dr(2)
DO ig = igmin, igmax
i = map(ig, 1)
x = (REAL(ig, dp) + REAL(cubecenter(1), dp))*dr(1)
point(1, 1) = x; point(2, 1) = y; point(3, 1) = z
point(1, 2) = x; point(2, 2) = y2; point(3, 2) = z
point(1, 3) = x; point(2, 3) = y; point(3, 3) = z2
point(1, 4) = x; point(2, 4) = y2; point(3, 4) = z2
!
res(1) = (point(ir, 1) - rb(ir)) - gauge(i, j, k)
res(2) = (point(ir, 2) - rb(ir)) - gauge(i, j2, k)
res(3) = (point(ir, 3) - rb(ir)) - gauge(i, j, k2)
res(4) = (point(ir, 4) - rb(ir)) - gauge(i, j2, k2)
!
grid_buf(i, j, k) = grid_buf(i, j, k) + grid(i, j, k)*res(1)
grid_buf(i, j2, k) = grid_buf(i, j2, k) + grid(i, j2, k)*res(2)
grid_buf(i, j, k2) = grid_buf(i, j, k2) + grid(i, j, k2)*res(3)
grid_buf(i, j2, k2) = grid_buf(i, j2, k2) + grid(i, j2, k2)*res(4)
ENDDO
ENDDO
ENDDO
END SUBROUTINE collocate_gauge_ortho

! **************************************************************************************************
!> \brief ...
!> \param current_env ...
Expand Down

0 comments on commit e4387f1

Please sign in to comment.