Skip to content

Commit

Permalink
Refactor RPA (#272)
Browse files Browse the repository at this point in the history 
* Refactor Laplace-MP2

* Change input structure for cubic scaling.

* Reactivate old IM_TIME RPA inputs

* Debugging Input, Laplace MP2, Remove unused Var. in mp2_laplace

* Refactor rpa_ri_gpw.F

* Add intent and pure in rpa_ri_gpw.F

* Refactor GW in rpa_ri_gpw.F

* Split/ Debug rpa_ri_gpw.F

* Use Fortran intrinsics in RPA/ Laplace-MP2

* Make routines PURE

* Do further factorization and cleaning up

* Make pretty

* Remove ri_metric_gw

* Move axk routines

* Cleanup

* Debug rpa_axk

* Remove RI_G0W0%STOP_CRIT

* Remove RI_G0W0%SCALING

* Adjust tolerances of some regtests

* Remove RI_G0W0%MIX_EXCHANGE

* Remove RI_G0W0%ATOMS

* Remove Contour deformation (It was buggy)

* Remove WFC_GPW%DO_CHOLESKY_SUBGROUPS

* Remove MULTIPOLE_TWO_CENT_INT

* Remove MINIMAL_GAP

* Set default group size for integral calc to 1

* Make pretty, return to ol do_regtest

* Push size/end/start array into calculate_BIb_C_2D

* Remove unused and move later used variables

* Fix line truncation

* Make pretty

* Debug rpa_gw and mp2_ri_gpw

* Remove unused argument

* Copy GW routines back to rpa_ri_gpw.F
  • Loading branch information
@fstein93 @oschuett
fstein93 authored and oschuett committed Apr 4, 2019
1 parent 48dde32 commit 7ea2c6e
Show file tree
Hide file tree
Showing 25 changed files with 7,932 additions and 4,444 deletions.
244 changes: 51 additions & 193 deletions src/input_cp2k_mp2.F
View file

Large diffs are not rendered by default.

16 changes: 4 additions & 12 deletions src/mp2.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1263,8 +1263,8 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
myfun, myfun_aux, myfun_prim, &
n_rep_hf, nmo, ns
LOGICAL :: charge_constrain_tmp, do_admm_rpa, &
do_hfx, do_mix_exchange, do_ri_Sigma_x
REAL(KIND=dp) :: eh1, ehfx, frac_exx, hfx_fraction, t1, t2
do_hfx, do_ri_Sigma_x
REAL(KIND=dp) :: eh1, ehfx, hfx_fraction, t1, t2
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_Sigma_x_minus_vxc_gw
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_type), POINTER :: mo_coeff
Expand Down Expand Up @@ -1444,12 +1444,10 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
do_mix_exchange = mp2_env%ri_g0w0%mix_exchange
! in most cases, we calculate the exchange self-energy here. But if we do only RI for
! the exchange self-energy, we do not calculate exchange here
ehfx = 0.0_dp
IF (.NOT. do_ri_Sigma_x .OR. do_mix_exchange) THEN
IF (.NOT. do_ri_Sigma_x) THEN
! add here HFX (=Sigma_exchange) to matrix_sigma_x_minus_vxc
DO irep = 1, n_rep_hf
Expand All @@ -1476,14 +1474,8 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
CALL admm_mo_merge_ks_matrix(qs_env)
END IF
IF (do_mix_exchange) THEN
frac_exx = mp2_env%ri_g0w0%frac_exx
ELSE
frac_exx = 1.0_dp
END IF
DO ispin = 1, dft_control%nspins
CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, frac_exx)
CALL dbcsr_add(matrix_sigma_x_minus_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
END DO
CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, mo_coeff_b)
Expand Down
46 changes: 3 additions & 43 deletions src/mp2_gpw.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
dimen_RI, group_size_3c, gw_corr_lev_occ, gw_corr_lev_occ_beta, gw_corr_lev_virt, &
gw_corr_lev_virt_beta, handle, homo, homo_beta, i, i_multigrid, local_unit_nr, &
my_group_L_end, my_group_L_size, my_group_L_start, n_multigrid, nelectron, &
nelectron_beta, nmo, nspins, ri_metric, ri_metric_gw
nelectron_beta, nmo, nspins, ri_metric
INTEGER, ALLOCATABLE, DIMENSION(:) :: ends_array, ends_B_all, ends_B_occ_bse, &
ends_B_virt_bse, ends_B_virtual, ends_B_virtual_beta, sizes_array, sizes_B_all, &
sizes_B_occ_bse, sizes_B_virt_bse, sizes_B_virtual, sizes_B_virtual_beta, starts_array, &
Expand Down Expand Up @@ -355,7 +355,6 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T

! in imaginary time, we do overlap metric
mp2_env%ri_metric = ri_overlap
mp2_env%ri_g0w0%ri_metric = ri_overlap

END IF

Expand Down Expand Up @@ -402,9 +401,6 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
gw_corr_lev_virt_beta = mp2_env%ri_g0w0%corr_mos_virt_beta
END IF

! metric for GW matrix elements B^nm_P
ri_metric_gw = mp2_env%ri_g0w0%ri_metric

! and the array of mos
ALLOCATE (Eigenval(dimen))
Eigenval(:) = mo_eigenvalues(:)
Expand All @@ -413,13 +409,6 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
Eigenval_beta(:) = mo_eigenvalues_beta(:)
ENDIF

IF (mp2_env%minimal_gap > 0.0_dp) THEN
CALL shift_eigenvalues(Eigenval, mp2_env%minimal_gap, homo, dimen)
IF (nspins == 2) THEN
CALL shift_eigenvalues(Eigenval_beta, mp2_env%minimal_gap, homo_beta, dimen)
END IF
END IF

! for imag. time, we do not need this
IF (.NOT. do_im_time) THEN

Expand Down Expand Up @@ -515,7 +504,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
do_bse, starts_B_all, sizes_B_all, ends_B_all, gw_corr_lev_occ, gw_corr_lev_virt, &
do_im_time, do_mao, do_kpoints_cubic_RPA, &
mat_3c_overl_int, mat_3c_overl_int_mao_for_occ, mat_3c_overl_int_mao_for_virt, &
mao_coeff_occ, mao_coeff_virt, ri_metric, ri_metric_gw, &
mao_coeff_occ, mao_coeff_virt, ri_metric, &
starts_B_occ_bse, sizes_B_occ_bse, ends_B_occ_bse, &
starts_B_virt_bse, sizes_B_virt_bse, ends_B_virt_bse, &
BIb_C_beta, BIb_C_gw_beta, ends_B_virtual_beta, sizes_B_virtual_beta, starts_B_virtual_beta, &
Expand All @@ -538,7 +527,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
do_im_time, do_mao, do_kpoints_cubic_RPA, &
mat_3c_overl_int, mat_3c_overl_int_mao_for_occ, &
mat_3c_overl_int_mao_for_virt, mao_coeff_occ, mao_coeff_virt, &
ri_metric, ri_metric_gw, &
ri_metric, &
starts_B_occ_bse, sizes_B_occ_bse, ends_B_occ_bse, &
starts_B_virt_bse, sizes_B_virt_bse, ends_B_virt_bse)
Expand Down Expand Up @@ -3837,33 +3826,4 @@ SUBROUTINE get_eps_old(dft_control, eps_pgf_orb_old, eps_rho_rspace_old, eps_gvg

END SUBROUTINE

! **************************************************************************************************
!> \brief ...
!> \param Eigenval ...
!> \param minimal_gap ...
!> \param homo ...
!> \param dimen ...
! **************************************************************************************************
SUBROUTINE shift_eigenvalues(Eigenval, minimal_gap, homo, dimen)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval
REAL(KIND=dp) :: minimal_gap
INTEGER :: homo, dimen

INTEGER :: n_level
REAL(KIND=dp) :: gap, minimal_occ, minimal_virt

gap = Eigenval(homo+1)-Eigenval(homo)
IF (gap < minimal_gap) THEN
minimal_occ = 0.5_dp*(Eigenval(homo)+Eigenval(homo+1)-minimal_gap)
minimal_virt = 0.5_dp*(Eigenval(homo)+Eigenval(homo+1)+minimal_gap)
DO n_level = 1, homo
IF (Eigenval(n_level) > minimal_occ) Eigenval(n_level) = minimal_occ
END DO
DO n_level = homo+1, dimen
IF (Eigenval(n_level) < minimal_virt) Eigenval(n_level) = minimal_virt
END DO
END IF

END SUBROUTINE

END MODULE mp2_gpw
51 changes: 16 additions & 35 deletions src/mp2_laplace.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,9 +82,8 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v
CHARACTER(LEN=*), PARAMETER :: routineN = 'laplace_minimax_approx', &
routineP = moduleN//':'//routineN

INTEGER :: avirt, handle, i_global, iiB, iocc, j_global, jjB, jquad, my_num_dgemm_call, &
ncol_local, nrow_local, number_of_rec, number_of_rec_beta, number_of_send, &
number_of_send_beta
INTEGER :: avirt, handle, i_global, iiB, iocc, jjB, jquad, my_num_dgemm_call, ncol_local, &
nrow_local, number_of_rec, number_of_rec_beta, number_of_send, number_of_send_beta
INTEGER, ALLOCATABLE, DIMENSION(:) :: map_rec_size, map_rec_size_beta, &
map_send_size, map_send_size_beta, &
RPA_proc_map
Expand Down Expand Up @@ -129,14 +128,10 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v
CALL get_exp_minimax_coeff(num_integ_points, E_Range, awj)

ALLOCATE (aj(num_integ_points))
aj = 0.0_dp
ALLOCATE (wj(num_integ_points))
wj = 0.0_dp

DO jquad = 1, num_integ_points
aj(jquad) = awj(jquad)
wj(jquad) = awj(jquad+num_integ_points)
END DO
aj(:) = awj(1:num_integ_points)
wj(:) = awj(num_integ_points+1:2*num_integ_points)

DEALLOCATE (awj)

Expand Down Expand Up @@ -180,21 +175,18 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v
!XXX CALL cp_fm_get_info(matrix=fm_mat_G,&
!XXX nrow_local=nrow_local,&
!XXX ncol_local=ncol_local,&
!XXX row_indices=row_indices,&
!XXX col_indices=col_indices)
!XXX row_indices=row_indices)
!XXX

! get info of fm_mat_S
CALL cp_fm_get_info(matrix=fm_mat_S, &
nrow_local=nrow_local, &
ncol_local=ncol_local, &
row_indices=row_indices, &
col_indices=col_indices)
row_indices=row_indices)

! update G matrix with the new value of w and a
IF (first_cycle) THEN
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)

Expand All @@ -210,7 +202,6 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v
END DO
ELSE
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)

Expand Down Expand Up @@ -250,20 +241,17 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v
!XXX CALL cp_fm_get_info(matrix=fm_mat_G_beta,&
!XXX nrow_local=nrow_local,&
!XXX ncol_local=ncol_local,&
!XXX row_indices=row_indices,&
!XXX col_indices=col_indices)
!XXX row_indices=row_indices)
!XXX

! the same for the beta spin
CALL cp_fm_get_info(matrix=fm_mat_S_beta, &
nrow_local=nrow_local, &
ncol_local=ncol_local, &
row_indices=row_indices, &
col_indices=col_indices)
row_indices=row_indices)
! update G matrix with the new value of w and a
IF (first_cycle) THEN
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)

Expand All @@ -280,7 +268,6 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v
END DO
ELSE
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)

Expand Down Expand Up @@ -315,26 +302,20 @@ SUBROUTINE laplace_minimax_approx(Emp2, para_env, para_env_RPA, unit_nr, homo, v

END IF

! get info of fm_mat_Q
CALL cp_fm_get_info(matrix=fm_mat_Q, &
nrow_local=nrow_local, &
ncol_local=ncol_local, &
row_indices=row_indices, &
col_indices=col_indices)

! calcualte the trace of the product Q*Q
! calculate the trace of the product Q*Q
trace_XX = 0.0_dp
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)
IF (my_open_shell) THEN
trace_XX = trace_XX+fm_mat_Q%local_data(iiB, jjB)*fm_mat_Q_beta%local_data(iiB, jjB)
ELSE
trace_XX = trace_XX+fm_mat_Q%local_data(iiB, jjB)*fm_mat_Q%local_data(iiB, jjB)
END IF
IF (my_open_shell) THEN
DO jjB = 1, ncol_local
trace_XX = trace_XX+DOT_PRODUCT(fm_mat_Q%local_data(:, jjB), fm_mat_Q_beta%local_data(:, jjB))
END DO
END DO
ELSE
trace_XX = NORM2(fm_mat_Q%local_data)
trace_XX = trace_XX*trace_XX
END IF

Emp2 = Emp2-trace_XX

Expand Down

0 comments on commit 7ea2c6e

Please sign in to comment.