Skip to content

Commit

Permalink
MP2: Refactor potential calculation GPW
Browse files Browse the repository at this point in the history
  • Loading branch information
@dev-zero
Frederick Stein authored and dev-zero committed Jul 5, 2019
1 parent 3659329 commit 554b7db
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 49 deletions.
17 changes: 7 additions & 10 deletions src/mp2_cphf.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ MODULE mp2_cphf
USE machine, ONLY: m_walltime
USE mathconstants, ONLY: fourpi
USE message_passing, ONLY: mp_sum
USE mp2_ri_gpw_in, ONLY: calc_potential_gpw
USE mp2_types, ONLY: mp2_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
Expand Down Expand Up @@ -152,8 +153,7 @@ SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
LOGICAL :: alpha_beta, do_dynamic_load_balancing, &
do_hfx, hfx_treat_lsd_in_core, &
restore_p_screen, use_virial
REAL(KIND=dp) :: e_hartree, factor, out_alpha, &
pair_energy, tot_rho_r
REAL(KIND=dp) :: e_hartree, factor, out_alpha, tot_rho_r
REAL(KIND=dp), DIMENSION(3, 3) :: h_stress
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
Expand Down Expand Up @@ -776,11 +776,8 @@ SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
total_rho=tot_rho_r, &
ks_env=ks_env, &
soft_valid=.FALSE.)
! calculate the MP2 potential
CALL pw_transfer(rho_r%pw, rho_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw)
CALL pw_transfer(pot_g%pw, rho_r%pw)
CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol)

CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g)

! calculate core forces
CALL integrate_v_core_rspace(rho_r, qs_env)
Expand Down Expand Up @@ -832,7 +829,7 @@ SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
CALL pw_transfer(rho_r%pw, rho_g%pw)
! don't forget the core density
CALL pw_axpy(rho_core%pw, rho_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, vhartree=pot_g%pw)
! finally update virial with the volume contribution
e_hartree = pw_integral_ab(temp_pw_g%pw, pot_g%pw)
Expand Down Expand Up @@ -964,7 +961,7 @@ SUBROUTINE cphf_like_update(qs_env, para_env, homo, virtual, dimen, &
jjB, ncol_local, nrow_local, ns
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
LOGICAL :: alpha_beta, my_recalc_hfx_integrals
REAL(KIND=dp) :: ehfx, pair_energy, total_rho
REAL(KIND=dp) :: ehfx, total_rho
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_work_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d
TYPE(qs_ks_env_type), POINTER :: ks_env
Expand Down Expand Up @@ -1215,7 +1212,7 @@ SUBROUTINE cphf_like_update(qs_env, para_env, homo, virtual, dimen, &
total_rho=total_rho, &
ks_env=ks_env)
! and calculate potential
CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, vhartree=pot_g%pw)
CALL pw_transfer(pot_g%pw, rho_r%pw)
CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol)
! integrate the potential
Expand Down
113 changes: 78 additions & 35 deletions src/mp2_ri_gpw_in.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,8 @@ MODULE mp2_ri_gpw_in
USE input_constants, ONLY: do_eri_gpw,&
do_eri_mme,&
do_eri_os,&
do_hfx_potential_coulomb,&
do_hfx_potential_id,&
ri_coulomb,&
ri_overlap
USE kinds, ONLY: dp
Expand Down Expand Up @@ -154,7 +156,7 @@ MODULE mp2_ri_gpw_in

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_gpw_in'

PUBLIC :: mp2_ri_gpw_compute_in
PUBLIC :: mp2_ri_gpw_compute_in, calc_potential_gpw

CONTAINS

Expand Down Expand Up @@ -313,13 +315,13 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
INTEGER :: color_sub_3c, cut_memory, cut_RI, eri_method, gw_corr_lev_total, handle, handle2, &
handle3, handle4, i, i_counter, itmp(2), j, LLL, max_row_col_local, &
max_row_col_local_beta, max_row_col_local_gw, max_row_col_local_occ_bse, &
max_row_col_local_virt_bse, my_B_all_end, my_B_all_size, my_B_all_start, &
my_B_occ_bse_end, my_B_occ_bse_size, my_B_occ_bse_start, my_B_size, my_B_size_beta, &
my_B_virt_bse_end, my_B_virt_bse_size, my_B_virt_bse_start, my_B_virtual_end, &
my_B_virtual_end_beta, my_B_virtual_start, my_B_virtual_start_beta, my_group_L_end, &
my_group_L_size, my_group_L_start, natom
INTEGER :: ngroup, ngroup_im_time, ngroup_im_time_P, nproc_sub, num_diff_blk, &
num_small_eigen, sqrt_max_bsize, virtual, virtual_beta
max_row_col_local_virt_bse, metric_potential, my_B_all_end, my_B_all_size, &
my_B_all_start, my_B_occ_bse_end, my_B_occ_bse_size, my_B_occ_bse_start, my_B_size, &
my_B_size_beta, my_B_virt_bse_end, my_B_virt_bse_size, my_B_virt_bse_start, &
my_B_virtual_end, my_B_virtual_end_beta, my_B_virtual_start, my_B_virtual_start_beta, &
my_group_L_end, my_group_L_size
INTEGER :: my_group_L_start, natom, ngroup, ngroup_im_time, ngroup_im_time_P, nproc_sub, &
num_diff_blk, num_small_eigen, sqrt_max_bsize, virtual, virtual_beta
INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_blk_sizes_split, local_atoms_for_mao_basis, &
local_atoms_RI, my_group_L_ends_im_time, my_group_L_sizes_im_time, &
my_group_L_starts_im_time, row_blk_sizes_RI_t, row_blk_sizes_RI_t_split, row_dist_RI_t, &
Expand All @@ -331,8 +333,7 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
INTEGER, DIMENSION(:), POINTER :: ao_blk_sizes
LOGICAL :: do_alpha_beta, do_kpoints_from_Gamma, &
memory_info
REAL(KIND=dp) :: cond_num, mem_for_iaK, pair_energy, &
wfn_size
REAL(KIND=dp) :: cond_num, mem_for_iaK, wfn_size
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: my_Lrows
TYPE(cp_eri_mme_param), POINTER :: eri_param
TYPE(cp_fm_type), POINTER :: fm_BIb_bse_ab, fm_BIb_bse_ij, fm_BIb_gw, &
Expand Down Expand Up @@ -710,6 +711,15 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
ELSE IF (eri_method == do_eri_gpw) THEN
SELECT CASE (ri_metric)
CASE (ri_coulomb)
metric_potential = do_hfx_potential_coulomb
CASE (ri_overlap)
metric_potential = do_hfx_potential_id
CASE DEFAULT
CPABORT("NIY")
END SELECT
DO LLL = my_group_L_start, my_group_L_end
i_counter = i_counter+1
Expand All @@ -719,19 +729,9 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
basis_type="RI_AUX", &
external_vector=my_Lrows(:, LLL-my_group_L_start+1))
CALL timeset(routineN//"_pot", handle3)
rho_r%pw%cr3d = psi_L%pw%cr3d
! in case we do Coulomb metric for RI, we need the Coulomb operator, but for RI with the
! overlap metric, we do not need the Coulomb operator
IF (ri_metric == ri_coulomb) THEN
CALL pw_transfer(rho_r%pw, rho_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw)
CALL pw_transfer(pot_g%pw, rho_r%pw)
END IF
CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol)
CALL timestop(handle3)
CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, metric_potential)
! and finally (K|mu nu)
CALL timeset(routineN//"_int", handle3)
Expand Down Expand Up @@ -2138,17 +2138,17 @@ SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_en
INTEGER :: best_group_size, color_L, comm_exchange, comm_L, dir, group_size, handle, &
handle2, handle3, i, i_counter, i_global, iatom, igrid_level, iiB, ikind, ipgf, iproc, &
iset, j_global, jjB, lb(3), LLL, location(3), my_ri_metric, na1, na2, natom, ncoa, &
ncol_local, nkind, nrow_local, nseta, nsgf, offset, proc_receive, proc_receive_static, &
proc_send, proc_send_static, proc_shift, rec_L_end, rec_L_size, rec_L_start, sgfa, &
strat_group_size, sub_sub_color, tp(3), ub(3)
iset, j_global, jjB, lb(3), LLL, location(3), metric_potential, my_ri_metric, na1, na2, &
natom, ncoa, ncol_local, nkind, nrow_local, nseta, nsgf, offset, proc_receive, &
proc_receive_static, proc_send, proc_send_static, proc_shift, rec_L_end, rec_L_size, &
rec_L_start, sgfa, strat_group_size, sub_sub_color, tp(3), ub(3)
INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, proc_map
INTEGER, DIMENSION(:), POINTER :: col_indices, la_max, la_min, npgfa, &
nsgfa, row_indices
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: map_it_here, my_do_im_time, &
my_external_matrix_struc
REAL(KIND=dp) :: min_mem_for_QK, pair_energy, rab2
REAL(KIND=dp) :: min_mem_for_QK, rab2
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egen_L, wf_vector
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab, L_external_col, L_local_col
REAL(KIND=dp), DIMENSION(3) :: ra, rab
Expand All @@ -2174,7 +2174,18 @@ SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_en
! default: Coulomb metric
my_ri_metric = ri_coulomb
IF (PRESENT(ri_metric)) my_ri_metric = ri_metric
metric_potential = do_hfx_potential_coulomb
IF (PRESENT(ri_metric)) THEN
my_ri_metric = ri_metric
SELECT CASE (ri_metric)
CASE (ri_coulomb)
metric_potential = do_hfx_potential_coulomb
CASE (ri_overlap)
metric_potential = do_hfx_potential_id
CASE DEFAULT
CPABORT("NIY")
END SELECT
END IF
my_external_matrix_struc = .FALSE.
IF (PRESENT(external_matrix_struc)) THEN
Expand Down Expand Up @@ -2254,14 +2265,7 @@ SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_en
CALL timeset(routineN//"_pot_lm", handle3)
rho_r%pw%cr3d = psi_L%pw%cr3d
! in case we do the overlap metric, we have no potential and hence no FT has to be carried out
IF (my_ri_metric == ri_coulomb) THEN
CALL pw_transfer(rho_r%pw, rho_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw)
CALL pw_transfer(pot_g%pw, rho_r%pw)
END IF
CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol)
CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, metric_potential)
NULLIFY (rs_v)
NULLIFY (rs_descs)
Expand Down Expand Up @@ -5733,4 +5737,43 @@ SUBROUTINE setup_group_L_im_time(my_group_L_starts_im_time, my_group_L_ends_im_t
END SUBROUTINE setup_group_L_im_time
! **************************************************************************************************
!> \brief ...
!> \param rho_r ...
!> \param rho_g ...
!> \param poisson_env ...
!> \param pot_g ...
!> \param potential_type ...
! **************************************************************************************************
SUBROUTINE calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, potential_type)
TYPE(pw_p_type), INTENT(IN) :: rho_r, rho_g
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_p_type), INTENT(IN) :: pot_g
INTEGER, INTENT(IN), OPTIONAL :: potential_type
CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_potential_gpw', &
routineP = moduleN//':'//routineN
INTEGER :: handle, my_potential_type
CALL timeset(routineN, handle)
my_potential_type = do_hfx_potential_coulomb
IF (PRESENT(potential_type)) THEN
my_potential_type = potential_type
END IF
! in case we do Coulomb metric for RI, we need the Coulomb operator, but for RI with the
! overlap metric, we do not need the Coulomb operator
IF (my_potential_type == do_hfx_potential_coulomb) THEN
CALL pw_transfer(rho_r%pw, rho_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, vhartree=pot_g%pw)
CALL pw_transfer(pot_g%pw, rho_r%pw)
END IF
CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol)
CALL timestop(handle)
END SUBROUTINE calc_potential_gpw
END MODULE mp2_ri_gpw_in
7 changes: 3 additions & 4 deletions src/mp2_ri_grad.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ MODULE mp2_ri_grad
mp2_eri_3c_integrate,&
mp2_eri_deallocate_forces,&
mp2_eri_force
USE mp2_ri_gpw_in, ONLY: calc_potential_gpw
USE mp2_types, ONLY: integ_mat_buffer_type,&
integ_mat_buffer_type_2D,&
mp2_type
Expand Down Expand Up @@ -462,10 +463,8 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, dft_contr
external_vector=wf_vector)
rho_r%pw%cr3d = 0.50_dp*(rho_r%pw%cr3d+psi_L_beta%pw%cr3d)
ENDIF
CALL pw_transfer(rho_r%pw, rho_g%pw)
CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw)
CALL pw_transfer(pot_g%pw, rho_r%pw)
CALL pw_scale(rho_r%pw, rho_r%pw%pw_grid%dvol)

CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g)
CALL timestop(handle3)

IF (use_virial) THEN
Expand Down

0 comments on commit 554b7db

Please sign in to comment.