comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 06:05:29 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

GW/BSE: Fix indices in case of removed null-spaces for SVDs

Browse source
This commit is contained in:
Max Graml 2025-03-04 15:55:46 +01:00 • committed by GitHub
parent cb33178bc1
commit 0d4bf46ef2
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
7 changed files with 163 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

12
src/bse_print.F
View file

@ -765,8 +765,10 @@ CONTAINS
DO i = 1, nrow_local
DO j = 1, ncol_local
IF (ABS(fm%local_data(i, j)) > thresh) THEN
WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
ABS(fm%local_data(i, j))
IF (unit_nr > 0) THEN
WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
ABS(fm%local_data(i, j))
END IF
END IF
END DO
END DO
@ -774,8 +776,10 @@ CONTAINS
DO i = 1, nrow_local
DO j = 1, ncol_local
IF (ABS(fm%local_data(i, j)) > thresh) THEN
WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
fm%local_data(i, j)
IF (unit_nr > 0) THEN
WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
fm%local_data(i, j)
END IF
END IF
END DO
END DO

29
src/bse_util.F
View file

@ -1726,10 +1726,10 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_multipoles_mo'
INTEGER :: handle, idir, n_multipole, n_occ, &
n_virt, nao
n_virt, nao, nmo_mp2
REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
TYPE(cp_fm_struct_type), POINTER :: fm_struct_mp_ab_trunc, fm_struct_mp_ai_trunc, &
fm_struct_mp_ij_trunc, fm_struct_multipoles_ao, fm_struct_nao_nao
fm_struct_mp_ij_trunc, fm_struct_multipoles_ao, fm_struct_nao_nmo, fm_struct_nmo_nmo
TYPE(cp_fm_type) :: fm_work
TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_multipole_per_dir
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_multipole, matrix_s
@ -1773,9 +1773,11 @@ CONTAINS
! Now we transform them to MO
! n_occ is the number of occupied MOs, nao the number of all AOs
! Writing homo to n_occ instead if nmo,
! takes care of ADDED_MOS, which would overwrite nmo, if invoked
! takes care of ADDED_MOS, which would overwrite nmo of qs_env-mos, if invoked
CALL get_mo_set(mo_set=mos(1), homo=n_occ, nao=nao)
n_virt = nao - n_occ
! Takes into account removed nullspace values from SVD
nmo_mp2 = mo_coeff(1)%matrix_struct%ncol_global
n_virt = nmo_mp2 - n_occ
! At the end, we need four different layouts of matrices in this multiplication, e.g. for a dipole:
! D_pq = full multipole matrix for occupied and unoccupied
@ -1787,9 +1789,12 @@ CONTAINS
! fm_work = <\mu|\vec{r}|\nu> C_{\nu q} EQ.II
! Struct for the full multipole matrix
CALL cp_fm_struct_create(fm_struct_nao_nao, &
CALL cp_fm_struct_create(fm_struct_nao_nmo, &
fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
nao, nao)
nao, nmo_mp2)
CALL cp_fm_struct_create(fm_struct_nmo_nmo, &
fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
nmo_mp2, nmo_mp2)
! At the very end, we copy the multipoles corresponding to truncated BSE indices in i and a
CALL cp_fm_struct_create(fm_struct_mp_ai_trunc, para_env_BSE, &
@ -1812,19 +1817,20 @@ CONTAINS
! Need another temporary matrix to store intermediate result from right multiplication
! D = C_{mu a} <\mu|\vec{r}|\nu> C_{\nu i}
CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nao, name="multipole_work")
CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nmo, name="multipole_work")
CALL cp_fm_set_all(fm_work, 0.0_dp)
DO idir = 1, n_multipole
! Create the full multipole matrix per direction
CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nao_nao, name="multipoles_mo")
CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nmo_nmo, name="multipoles_mo")
CALL cp_fm_set_all(fm_multipole_per_dir(idir), 0.0_dp)
! Fill final (MO) multipole matrix
CALL cp_dbcsr_sm_fm_multiply(matrix_multipole(idir)%matrix, mo_coeff(1), &
fm_work, ncol=nao)
fm_work, ncol=nmo_mp2)
! Now obtain the multipoles by the final multiplication;
! We do that inside the loop to obtain multipoles per axis for print
CALL parallel_gemm('T', 'N', nao, nao, nao, 1.0_dp, mo_coeff(1), fm_work, 0.0_dp, fm_multipole_per_dir(idir))
CALL parallel_gemm('T', 'N', nmo_mp2, nmo_mp2, nao, 1.0_dp, mo_coeff(1), fm_work, 0.0_dp, fm_multipole_per_dir(idir))
! Truncate full matrix to the BSE indices
! D_ai
CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
@ -1863,7 +1869,8 @@ CONTAINS
CALL cp_fm_struct_release(fm_struct_mp_ai_trunc)
CALL cp_fm_struct_release(fm_struct_mp_ij_trunc)
CALL cp_fm_struct_release(fm_struct_mp_ab_trunc)
CALL cp_fm_struct_release(fm_struct_nao_nao)
CALL cp_fm_struct_release(fm_struct_nao_nmo)
CALL cp_fm_struct_release(fm_struct_nmo_nmo)
DO idir = 1, n_multipole
CALL cp_fm_release(fm_multipole_per_dir(idir))
END DO

2
src/qs_scf_methods.F
View file

@ -391,6 +391,8 @@ CONTAINS
mo_coeff)
CALL cp_fm_release(work_red2)
ELSE
CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, matrix_ks_fm, ortho, 0.0_dp, work)
CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, matrix_ks_fm)
IF (do_level_shift) &
CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm)

34
src/rpa_gw_sigma_x.F
View file

@ -123,9 +123,9 @@ CONTAINS
CHARACTER(4) :: occ_virt
CHARACTER(LEN=40) :: line
INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, handle, homo, &
homo_reduced_bse, homo_startindex_bse, i_img, ikp, irep, ispin, iunit, myfun, myfun_aux, &
myfun_prim, n_level_gw, n_level_gw_ref, n_rep_hf, nkp, nkp_Sigma, nmo, nspins, print_exx, &
virtual_reduced_bse, virtual_startindex_bse
homo_reduced_bse, homo_startindex_bse, i_img, ikp, irep, ispin, iunit, max_corr_lev_occ, &
max_corr_lev_virt, myfun, myfun_aux, myfun_prim, n_level_gw, n_level_gw_ref, n_rep_hf, &
nkp, nkp_Sigma, nmo, nspins, print_exx, virtual_reduced_bse, virtual_startindex_bse
LOGICAL :: calc_ints, charge_constrain_tmp, do_admm_rpa, do_hfx, do_kpoints_cubic_RPA, &
do_kpoints_from_Gamma, do_ri_Sigma_x, really_read_line
REAL(KIND=dp) :: E_GAP_GW, E_HOMO_GW, E_LUMO_GW, eh1, ehfx, eigval_dft, eigval_hf_at_dft, &
@ -558,6 +558,12 @@ CONTAINS
gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
! If SVD is used to invert overlap matrix (for CHOLESKY OFF), some MOs are removed
! Therefore, setting the number of gw_corr_lev_virt simply to dimen - homo leads to index problems
! Instead, we take into account the removed MOs
max_corr_lev_occ = homo
max_corr_lev_virt = nmo - homo
! If BSE is invoked, manipulate corrected MO number
IF (mp2_env%bse%do_bse) THEN
! Logic: If cutoff is negative, all MOs are included in BSE, i.e. we need to correct them all
@ -582,7 +588,7 @@ CONTAINS
! Obtain indices from DFT if gw_corr... are set to -2
CALL determine_cutoff_indices(mo_eigenvalues, &
homo, dimen - homo, &
homo, max_corr_lev_virt, &
homo_reduced_bse, virtual_reduced_bse, &
homo_startindex_bse, virtual_startindex_bse, &
mp2_env)
@ -599,8 +605,8 @@ CONTAINS
! if requested number of occ/virt levels for correction either exceed the number of
! occ/virt levels or the requested number is negative, default to correct all
! occ/virt level energies
IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = homo
IF (gw_corr_lev_virt > dimen - homo .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = dimen - homo
IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = max_corr_lev_occ
IF (gw_corr_lev_virt > max_corr_lev_virt .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = max_corr_lev_virt
IF (ispin == 1) THEN
mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt
@ -842,7 +848,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_sigma_x_minus_vxc_to_MO_basis'
INTEGER :: dimen, gw_corr_lev_occ, gw_corr_lev_virt, handle, homo, i_global, iiB, ikp, &
ispin, j_global, jjB, ncol_local, nkp, nrow_local, nspins
ispin, j_global, jjB, max_corr_lev_occ, max_corr_lev_virt, ncol_local, nkp, nrow_local, &
nspins
INTEGER, DIMENSION(2) :: kp_range
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
REAL(KIND=dp) :: imval, reval
@ -942,17 +949,22 @@ CONTAINS
CALL para_env%sum(vec_Sigma_x_minus_vxc_gw)
CALL para_env%sum(vec_Sigma_x_minus_vxc_gw_im)
! also adjust in the case of kpoints too big gw_corr_lev_occ and gw_corr_lev_virt
DO ispin = 1, nspins
CALL get_mo_set(mo_set=kpoints%kp_env(1)%kpoint_env%mos(ispin, 1), &
homo=homo, nao=dimen)
! If SVD is used to invert overlap matrix (for CHOLESKY OFF), some MOs are removed
! Therefore, setting the number of gw_corr_lev_virt simply to dimen - homo leads to index problems
! Instead, we take into account the removed MOs
max_corr_lev_occ = homo
max_corr_lev_virt = nmo - homo
gw_corr_lev_occ = mp2_env%ri_g0w0%corr_mos_occ
gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
! if corrected occ/virt levels exceed the number of occ/virt levels,
! if corrected occ/virt levels exceed the number of occ/virt levels or are negative,
! correct all occ/virt level energies
IF (gw_corr_lev_occ > homo) gw_corr_lev_occ = homo
IF (gw_corr_lev_virt > dimen - homo) gw_corr_lev_virt = dimen - homo
IF (gw_corr_lev_occ > homo .OR. gw_corr_lev_occ < 0) gw_corr_lev_occ = max_corr_lev_occ
IF (gw_corr_lev_virt > max_corr_lev_virt .OR. gw_corr_lev_virt < 0) gw_corr_lev_virt = max_corr_lev_virt
IF (ispin == 1) THEN
mp2_env%ri_g0w0%corr_mos_occ = gw_corr_lev_occ
mp2_env%ri_g0w0%corr_mos_virt = gw_corr_lev_virt

34
tests/QS/regtest-bse/BASIS_DIFFUSE Normal file
View file

@ -0,0 +1,34 @@
H def2-SVP-custom-diffuse
6
1 0 0 1 1
1.9622572 0.13796524
1 0 0 1 1
0.44453796 0.47831935
1 0 0 1 1
0.3 1.0000000
1 0 0 1 1
1.0 0.1000000
1 0 0 1 1
0.95 0.101
1 1 1 1 1
0.8000000 1.0000000
H def2-SVP-RIFIT-custom
9
1 0 0 1 1
9.33521609 .64609379
1 0 0 1 1
1.86110704 1.37005241
1 0 0 1 1
.59512466 1.0000000
1 0 0 1 1
.26448099 1.0000000
1 1 1 1 1
2.45249821 .13284956
1 1 1 1 1
1.35403830 1.45276215
1 1 1 1 1
.59522394 1.0000000
1 2 2 1 1
1.58163766 1.49367397
1 2 2 1 1
.62743960 -.01070690

77
tests/QS/regtest-bse/TDA_H2_PBE_G0W0_CHOLESKY_OFF.inp Normal file
View file

@ -0,0 +1,77 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT BSE_H2_PBE_CHOLESKY_OFF
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END TIMINGS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_DIFFUSE
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 300
REL_CUTOFF 30
&END MGRID
&POISSON
PERIODIC NONE
PSOLVER MULTIPOLE
&END POISSON
&PRINT
&MO
ENERGIES T
&END MO
&END PRINT
&QS
METHOD GAPW
&END QS
&SCF
CHOLESKY OFF
EPS_SCF 1.0E-6
MAX_SCF 200
SCF_GUESS RESTART
&END SCF
&XC
&WF_CORRELATION
&RI_RPA
&GW
SELF_CONSISTENCY G0W0
&BSE
ENERGY_CUTOFF_EMPTY 60
NUM_PRINT_EXC -1
TDA ON
&END BSE
&END GW
&END RI_RPA
&END WF_CORRELATION
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6 6 6
PERIODIC NONE
&END CELL
&COORD
H 0 0.0 0.0
H 0 0.0 0.74144
&END COORD
&KIND H
BASIS_SET def2-SVP-custom-diffuse
BASIS_SET RI_AUX def2-SVP-RIFIT-custom
POTENTIAL ALL
&END KIND
&PRINT
&ATOMIC_COORDINATES ON
&END ATOMIC_COORDINATES
&END PRINT
&TOPOLOGY
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

1
tests/QS/regtest-bse/TEST_FILES.toml
View file

@ -19,6 +19,7 @@
"BSE_H2O_PBE_evGW0_custom_screening.inp" = [{matcher="M110", tol=2e-05, ref=14.9538}]
"BSE_H2O_PBE_evGW0_spectra.inp" = [{matcher="M124", tol=2e-02, ref=25.3518},
{matcher="M125", tol=2e-02, ref=0.7006}]
"TDA_H2_PBE_G0W0_CHOLESKY_OFF.inp" = [{matcher="M109", tol=2e-05 , ref=16.9750}]
#EOF
#Low convergence criteria for evGW make a higher threshold necessary for these tests
#Logic: Absolute errors (for G0W0/evGW0) should be within 1e-4, i.e. relative errors are