Skip to content

Commit

Permalink
grid: Move radius calculation out of integrate_pgf_product
Browse files Browse the repository at this point in the history
  • Loading branch information
@oschuett
oschuett committed Feb 14, 2020
1 parent 5dc668a commit e1bf4e7
Show file tree
Hide file tree
Showing 8 changed files with 146 additions and 124 deletions.
61 changes: 9 additions & 52 deletions src/grid/grid_integrate.F
View file
Expand Up @@ -3,7 +3,6 @@
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
MODULE grid_integrate
USE ao_util, ONLY: exp_radius_very_extended
USE cell_types, ONLY: cell_type
USE cube_utils, ONLY: compute_cube_center,&
cube_info_type,&
Expand Down Expand Up @@ -71,16 +70,13 @@ END SUBROUTINE call_to_xyz_to_vab
!> \param pab ...
!> \param o1 ...
!> \param o2 ...
!> \param eps_gvg_rspace ...
!> \param radius ...
!> \param calculate_forces ...
!> \param hdab ...
!> \param hadb ...
!> \param force_a ...
!> \param force_b ...
!> \param compute_tau ...
!> \param map_consistent ...
!> \param collocate_rho0 ...
!> \param rpgf0_s ...
!> \param use_virial ...
!> \param my_virial_a ...
!> \param my_virial_b ...
Expand All @@ -92,10 +88,10 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, &
lb_max, zetb, lb_min, &
ra, rab, rsgrid, cell, &
cube_info, hab, pab, o1, o2, &
eps_gvg_rspace, &
radius, &
calculate_forces, hdab, hadb, force_a, force_b, &
compute_tau, map_consistent, &
collocate_rho0, rpgf0_s, use_virial, my_virial_a, &
compute_tau, &
use_virial, my_virial_a, &
my_virial_b, a_hdab, use_subpatch, subpatch_pattern)
INTEGER, INTENT(IN) :: la_max
Expand All @@ -111,15 +107,13 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, &
REAL(KIND=dp), DIMENSION(:, :), &
OPTIONAL, POINTER :: pab
INTEGER, INTENT(IN) :: o1, o2
REAL(KIND=dp), INTENT(IN) :: eps_gvg_rspace
REAL(KIND=dp), INTENT(IN) :: radius
LOGICAL, INTENT(IN) :: calculate_forces
REAL(KIND=dp), DIMENSION(:, :, :), &
OPTIONAL, POINTER :: hdab, hadb
REAL(KIND=dp), DIMENSION(3), &
INTENT(INOUT), OPTIONAL :: force_a, force_b
LOGICAL, INTENT(IN), OPTIONAL :: compute_tau, map_consistent, &
collocate_rho0
REAL(dp), INTENT(IN), OPTIONAL :: rpgf0_s
LOGICAL, INTENT(IN), OPTIONAL :: compute_tau
LOGICAL, INTENT(IN), OPTIONAL :: use_virial
REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: my_virial_a, my_virial_b
REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: a_hdab
Expand All @@ -136,15 +130,11 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, &
INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, &
ub_cube
INTEGER, DIMENSION(:), POINTER :: sphere_bounds
LOGICAL :: my_collocate_rho0, &
my_compute_tau, &
my_map_consistent, &
LOGICAL :: my_compute_tau, &
my_use_virial, &
subpatch_integrate
REAL(KIND=dp) :: axpm0, cutoff, &
der_a(3), der_b(3), f, ftza, ftzb, pabval, pg, &
prefactor, radius, rpg, zetp
REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp
REAL(KIND=dp) :: axpm0, der_a(3), der_b(3), f, ftza, ftzb, pabval, pg, prefactor, rpg, zetp
REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, roffset, rp
REAL(KIND=dp), DIMENSION(:, :, :), &
POINTER :: grid
Expand Down Expand Up @@ -183,19 +173,6 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, &
my_compute_tau = .FALSE.
ENDIF
! use identical radii for integrate and collocate ?
IF (PRESENT(map_consistent)) THEN
my_map_consistent = map_consistent
ELSE
my_map_consistent = .FALSE.
ENDIF
IF (PRESENT(collocate_rho0) .AND. PRESENT(rpgf0_s)) THEN
my_collocate_rho0 = collocate_rho0
ELSE
my_collocate_rho0 = .FALSE.
END IF
IF (calculate_forces) THEN
la_max_local = la_max + 1 ! needed for the derivative of the gaussian, unimportant which one
la_min_local = MAX(la_min - 1, 0) ! just in case the la_min,lb_min is not zero
Expand Down Expand Up @@ -225,29 +202,9 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, &
prefactor = EXP(-zeta*f*DOT_PRODUCT(rab, rab))
! *** position of the gaussian product
rap(:) = f*rab(:)
rbp(:) = rap(:) - rab(:)
rp(:) = ra(:) + rap(:) ! this is the gaussian center in real coordinates
rb(:) = ra(:) + rab(:)
IF (my_map_consistent) THEN ! still assumes that eps_gvg_rspace=eps_rho_rspace
cutoff = 1.0_dp
radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, ra=ra, rb=rb, rp=rp, &
zetp=zetp, eps=eps_gvg_rspace, prefactor=prefactor, cutoff=cutoff)
ELSE IF (my_collocate_rho0) THEN
cutoff = 0.0_dp
prefactor = 1.0_dp
radius = rpgf0_s
ELSE
cutoff = 1.0_dp
IF (PRESENT(pab)) THEN
radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, &
zetp, eps_gvg_rspace, prefactor, cutoff)
ELSE
radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, ra=ra, rb=rb, rp=rp, &
zetp=zetp, eps=eps_gvg_rspace, prefactor=prefactor, cutoff=cutoff)
ENDIF
ENDIF
IF (radius == 0.0_dp) THEN
RETURN
ENDIF
Expand Down
12 changes: 9 additions & 3 deletions src/hirshfeld_methods.F
View file
Expand Up @@ -543,7 +543,7 @@ SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv)
INTEGER, ALLOCATABLE, DIMENSION(:) :: cores
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: do_force
REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace
REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace, radius
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
Expand Down Expand Up @@ -610,11 +610,17 @@ SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv)
hab(1, 1) = 0.0_dp
force_a(:) = 0.0_dp
force_b(:) = 0.0_dp
!

radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
ra=ra, rb=ra, rp=ra, &
zetp=alpha, eps=eps_rho_rspace, &
pab=pab, o1=0, o2=0, & ! without map_consistent
prefactor=1.0_dp, cutoff=1.0_dp)

CALL integrate_pgf_product(0, alpha, 0, &
0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
rs_v, cell, pw_env%cube_info(1), hab, pab=pab, o1=0, o2=0, &
eps_gvg_rspace=eps_rho_rspace, calculate_forces=do_force, &
radius=radius, calculate_forces=do_force, &
force_a=force_a, force_b=force_b, use_virial=.FALSE., &
use_subpatch=.TRUE., subpatch_pattern=0_int_8)
fval(atom_a) = fval(atom_a) + hab(1, 1)*dvol*coef
Expand Down
18 changes: 11 additions & 7 deletions src/mp2_eri_gpw.F
View file
Expand Up @@ -9,6 +9,7 @@
!> 07.2019 separated from mp2_integrals.F [Frederick Stein]
! **************************************************************************************************
MODULE mp2_eri_gpw
USE ao_util, ONLY: exp_radius_very_extended
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_types, ONLY: gto_basis_set_type
USE cell_types, ONLY: cell_type,&
Expand Down Expand Up @@ -182,10 +183,10 @@ SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: map_it_here
REAL(KIND=dp) :: cutoff_old, relative_cutoff_old
REAL(KIND=dp) :: cutoff_old, radius, relative_cutoff_old
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab
REAL(KIND=dp), DIMENSION(3) :: ra, rab
REAL(KIND=dp), DIMENSION(3) :: ra
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, rpgfa, sphi_a, zeta
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
Expand Down Expand Up @@ -260,7 +261,6 @@ SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my
zeta => basis_set_a%zet

ra(:) = pbc(particle_set(iatom)%r, cell)
rab = 0.0_dp

DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
Expand Down Expand Up @@ -303,23 +303,27 @@ SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my
IF (map_it_here) THEN
DO ipgf = 1, npgfa(iset)
sgfa = first_sgfa(1, iset)

na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
na2 = ipgf*ncoset(la_max(iset))
igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))

radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
zetp=zeta(ipgf, iset), &
eps=dft_control%qs_control%eps_gvg_rspace, &
prefactor=1.0_dp, cutoff=1.0_dp)

CALL integrate_pgf_product( &
la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
ra=ra, rab=rab, &
ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
rsgrid=rs_v(igrid_level)%rs_grid, &
cell=cell, &
cube_info=pw_env_sub%cube_info(igrid_level), &
hab=I_tmp2, &
o1=na1 - 1, &
o2=0, &
map_consistent=.TRUE., &
eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, &
radius=radius, &
calculate_forces=.FALSE.)
END DO

Expand Down
29 changes: 19 additions & 10 deletions src/mp2_ri_grad.F
View file
Expand Up @@ -9,6 +9,7 @@
!> 10.2013 created [Mauro Del Ben]
! **************************************************************************************************
MODULE mp2_ri_grad
USE ao_util, ONLY: exp_radius_very_extended
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE basis_set_types, ONLY: gto_basis_set_type
Expand Down Expand Up @@ -171,11 +172,11 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
LOGICAL :: alpha_beta, map_it_here, skip_shell, &
use_virial
REAL(KIND=dp) :: cutoff_old, e_hartree, eps_filter, &
factor, omega, pair_energy, &
factor, omega, pair_energy, radius, &
relative_cutoff_old, total_rho
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: G_PQ_local, I_ab
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, my_virial_a, my_virial_b
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, pab, rpgfa, sphi_a, zeta
Expand Down Expand Up @@ -507,7 +508,6 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
zeta => basis_set_a%zet

ra(:) = pbc(particle_set(iatom)%r, cell)
rab = 0.0_dp

force_a(:) = 0.0_dp
force_b(:) = 0.0_dp
Expand Down Expand Up @@ -571,18 +571,23 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
na2 = ipgf*ncoset(la_max(iset))
igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))

radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
zetp=zeta(ipgf, iset), &
eps=dft_control%qs_control%eps_gvg_rspace, &
prefactor=1.0_dp, cutoff=1.0_dp)

CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
ra=ra, rab=rab, &
ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
rsgrid=rs_v(igrid_level)%rs_grid, &
cell=cell, &
cube_info=pw_env_sub%cube_info(igrid_level), &
hab=I_tmp2, &
pab=pab, &
o1=na1 - 1, &
o2=0, &
map_consistent=.TRUE., &
eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, &
radius=radius, &
calculate_forces=.TRUE., &
force_a=force_a, force_b=force_b, &
use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
Expand Down Expand Up @@ -744,7 +749,6 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
zeta => basis_set_a%zet

ra(:) = pbc(particle_set(iatom)%r, cell)
rab = 0.0_dp

force_a(:) = 0.0_dp
force_b(:) = 0.0_dp
Expand Down Expand Up @@ -821,18 +825,23 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
na2 = ipgf*ncoset(la_max(iset))
igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset))

radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
zetp=zeta(ipgf, iset), &
eps=dft_control%qs_control%eps_gvg_rspace, &
prefactor=1.0_dp, cutoff=1.0_dp)

CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), &
lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, &
ra=ra, rab=rab, &
ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
rsgrid=rs_v(igrid_level)%rs_grid, &
cell=cell, &
cube_info=pw_env_sub%cube_info(igrid_level), &
hab=I_tmp2, &
pab=pab, &
o1=na1 - 1, &
o2=0, &
map_consistent=.TRUE., &
eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, &
radius=radius, &
calculate_forces=.TRUE., &
force_a=force_a, force_b=force_b, &
use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
Expand Down

0 comments on commit e1bf4e7

Please sign in to comment.