Skip to content

Commit

Permalink
MP2/RPA: Added SVD as alternative for Cholesky
Browse files Browse the repository at this point in the history 
This is a preparation for the implementation of longrange potentials.
  • Loading branch information
@oschuett
Frederick Stein authored and oschuett committed Sep 18, 2019
1 parent 5566d17 commit 0a87912
Show file tree
Hide file tree
Showing 37 changed files with 1,320 additions and 122 deletions.
12 changes: 7 additions & 5 deletions src/fm/cp_fm_basic_linalg.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1127,8 +1127,7 @@ END SUBROUTINE cp_fm_scale
!> \param matrix ...
!> \param matrixt ...
!> \note
!> all matrix elements are transpose (see cp_fm_upper_to_half to symmetrise a matrix)
!> all matrix elements are transpose (see cp_fm_upper_to_half to symmetrize a matrix)
!> all matrix elements are transposed (see cp_fm_upper_to_half to symmetrise a matrix)
! **************************************************************************************************
SUBROUTINE cp_fm_transpose(matrix, matrixt)
TYPE(cp_fm_type), POINTER :: matrix, matrixt
Expand All @@ -1137,7 +1136,7 @@ SUBROUTINE cp_fm_transpose(matrix, matrixt)
routineP = moduleN//':'//routineN

INTEGER :: handle, ncol_global, &
nrow_global
nrow_global, ncol_globalt, nrow_globalt
REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, c
#if defined(__SCALAPACK)
INTEGER, DIMENSION(9) :: desca, descc
Expand All @@ -1149,7 +1148,10 @@ SUBROUTINE cp_fm_transpose(matrix, matrixt)
CPASSERT(ASSOCIATED(matrixt))
nrow_global = matrix%matrix_struct%nrow_global
ncol_global = matrix%matrix_struct%ncol_global
CPASSERT(nrow_global == ncol_global)
nrow_globalt = matrixt%matrix_struct%nrow_global
ncol_globalt = matrixt%matrix_struct%ncol_global
CPASSERT(nrow_global == ncol_globalt)
CPASSERT(nrow_globalt == ncol_global)

CALL timeset(routineN, handle)

Expand All @@ -1159,7 +1161,7 @@ SUBROUTINE cp_fm_transpose(matrix, matrixt)
#if defined(__SCALAPACK)
desca(:) = matrix%matrix_struct%descriptor(:)
descc(:) = matrixt%matrix_struct%descriptor(:)
CALL pdtran(nrow_global, ncol_global, 1.0_dp, a(1, 1), 1, 1, desca, 0.0_dp, c(1, 1), 1, 1, descc)
CALL pdtran(ncol_global, nrow_global, 1.0_dp, a(1, 1), 1, 1, desca, 0.0_dp, c(1, 1), 1, 1, descc)
#else
DO j = 1, ncol_global
DO i = 1, nrow_global
Expand Down
21 changes: 19 additions & 2 deletions src/input_cp2k_mp2.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,8 +219,7 @@ SUBROUTINE create_mp2_section(section)
CALL section_release(print_key)

CALL keyword_create( &
keyword, __LOCATION__, &
name="IM_TIME", &
keyword, __LOCATION__, name="IM_TIME", &
variants=(/"IMAG_TIME"/), &
description="Turns on cubic scaling RI-RPA, GW and Laplace-SOS-MP2 method using the imaginary time formalism. "// &
"If no IM_TIME section is present, the default parameters are used.", &
Expand All @@ -230,6 +229,24 @@ SUBROUTINE create_mp2_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DO_SVD", &
description="Wether to perform a singular value decomposition instead of the Cholesky decomposition "// &
"of the potential operator in the RI basis. Computationally expensive but numerically more stable. "// &
"It reduces the coomputational cost of some subsequent steps. Recommended when a longrange Coulomb "// &
"potential is employed for a much emproved numerical stability.", &
usage="DO_SVD .TRUE.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="EPS_SVD", &
description="Include all singular vectors with a singular value smaller than EPS_SVD. "// &
"Is used to provide info on smallest eigenvalues.", &
usage="EPS_SVD 1.0e-5", &
default_r_val=1.0E-3_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL create_mp2_direct(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
Expand Down
24 changes: 12 additions & 12 deletions src/mp2_gpw.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,10 +186,10 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_gpw_main', routineP = moduleN//':'//routineN

INTEGER :: blacs_grid_layout, color_sub, color_sub_3c, comm_sub, comm_sub_3c, dimen, &
dimen_RI, group_size_3c, group_size_P, 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
dimen_RI, dimen_RI_red, group_size_3c, group_size_P, 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
INTEGER, ALLOCATABLE, DIMENSION(:) :: ends_array_mc_t, starts_array_mc_t
LOGICAL :: blacs_repeatable, do_bse, do_dbcsr_t, do_im_time, do_kpoints_cubic_RPA, do_mao, &
my_do_gw, my_do_ri_mp2, my_do_ri_rpa, my_do_ri_sos_laplace_mp2, &
Expand Down Expand Up @@ -515,7 +515,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
IF (nspins == 2) THEN
! open shell case (RI) here the (ia|K) integrals are computed for both the alpha and beta components
CALL mp2_ri_gpw_compute_in( &
BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, dimen_RI, qs_env, &
BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, dimen_RI, dimen_RI_red, qs_env, &
para_env, para_env_sub, color_sub, dft_control, cell, particle_set, &
atomic_kind_set, qs_kind_set, mo_coeff, fm_matrix_L_RI_metric, nmo, homo, rho_r, rho_g, pot_g, &
mat_munu, mat_munu_mao_occ_virt, mat_munu_mao_virt_occ, sab_orb_sub, sab_orb_all, &
Expand All @@ -534,7 +534,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
ELSE
! closed shell case (RI)
CALL mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, &
dimen_RI, qs_env, para_env, para_env_sub, &
dimen_RI, dimen_RI_red, qs_env, para_env, para_env_sub, &
color_sub, dft_control, cell, particle_set, &
atomic_kind_set, qs_kind_set, mo_coeff, fm_matrix_L_RI_metric, &
nmo, homo, rho_r, rho_g, pot_g, &
Expand Down Expand Up @@ -729,7 +729,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
para_env, para_env_sub_RPA, color_sub, &
gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
mo_coeff, fm_matrix_L_RI_metric, kpoints, &
Eigenval, nmo, homo, dimen_RI, gw_corr_lev_occ, gw_corr_lev_virt, &
Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, &
unit_nr, my_do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_mao, do_bse, matrix_s, &
mao_coeff_occ, mao_coeff_virt, mao_coeff_occ_A, mao_coeff_virt_A, &
mat_munu, mat_dm_occ_local, mat_dm_virt_local, &
Expand All @@ -751,7 +751,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
para_env, para_env_sub_RPA, color_sub, &
gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
mo_coeff, fm_matrix_L_RI_metric, kpoints, &
Eigenval, nmo, homo, dimen_RI, gw_corr_lev_occ, gw_corr_lev_virt, &
Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, &
unit_nr, my_do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_mao, do_bse, matrix_s, &
mao_coeff_occ, mao_coeff_virt, mao_coeff_occ_A, mao_coeff_virt_A, &
mat_munu, mat_dm_occ_local, mat_dm_virt_local, &
Expand Down Expand Up @@ -800,21 +800,21 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
CALL mp2_ri_gpw_compute_en( &
Emp2, Emp2_Cou, Emp2_EX, BIb_C, mp2_env, para_env, para_env_sub, color_sub, &
gd_array, gd_B_virtual, &
Eigenval, nmo, homo, dimen_RI, unit_nr, calc_forces, calc_ex, &
Eigenval, nmo, homo, dimen_RI_red, unit_nr, calc_forces, calc_ex, &
open_shell_SS=.TRUE.)

! beta-beta component
CALL mp2_ri_gpw_compute_en( &
Emp2_BB, Emp2_Cou_BB, Emp2_EX_BB, BIb_C_beta, mp2_env, para_env, para_env_sub, color_sub, &
gd_array, gd_B_virtual_beta, &
Eigenval_beta, nmo, homo_beta, dimen_RI, unit_nr, calc_forces, calc_ex, &
Eigenval_beta, nmo, homo_beta, dimen_RI_red, unit_nr, calc_forces, calc_ex, &
open_shell_SS=.TRUE.)

! alpha-beta case
CALL mp2_ri_gpw_compute_en( &
Emp2_d_AB, Emp2_AB, Emp2_d2_AB, BIb_C, mp2_env, para_env, para_env_sub, color_sub, &
gd_array, gd_B_virtual, &
Eigenval, nmo, homo, dimen_RI, unit_nr, calc_forces, .FALSE., &
Eigenval, nmo, homo, dimen_RI_red, unit_nr, calc_forces, .FALSE., &
.FALSE., BIb_C_beta, homo_beta, Eigenval_beta, &
gd_B_virtual_beta)

Expand All @@ -823,7 +823,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
CALL mp2_ri_gpw_compute_en( &
Emp2, Emp2_Cou, Emp2_EX, BIb_C, mp2_env, para_env, para_env_sub, color_sub, &
gd_array, gd_B_virtual, &
Eigenval, nmo, homo, dimen_RI, unit_nr, calc_forces, calc_ex)
Eigenval, nmo, homo, dimen_RI_red, unit_nr, calc_forces, calc_ex)
END IF
! if we need forces time to calculate the MP2 non-separable contribution
! and start coputing the Lagrangian
Expand Down

0 comments on commit 0a87912

Please sign in to comment.