diff --git a/src/bse_print.F b/src/bse_print.F index 821b46ad27..be734350b9 100644 --- a/src/bse_print.F +++ b/src/bse_print.F @@ -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 diff --git a/src/bse_util.F b/src/bse_util.F index 673705fd3a..4969c07dff 100644 --- a/src/bse_util.F +++ b/src/bse_util.F @@ -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 diff --git a/src/qs_scf_methods.F b/src/qs_scf_methods.F index 41337d4bdc..d9cb09bcab 100644 --- a/src/qs_scf_methods.F +++ b/src/qs_scf_methods.F @@ -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) diff --git a/src/rpa_gw_sigma_x.F b/src/rpa_gw_sigma_x.F index 799af858f6..504a3638fe 100644 --- a/src/rpa_gw_sigma_x.F +++ b/src/rpa_gw_sigma_x.F @@ -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 diff --git a/tests/QS/regtest-bse/BASIS_DIFFUSE b/tests/QS/regtest-bse/BASIS_DIFFUSE new file mode 100644 index 0000000000..a3dfe78c2f --- /dev/null +++ b/tests/QS/regtest-bse/BASIS_DIFFUSE @@ -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 \ No newline at end of file diff --git a/tests/QS/regtest-bse/TDA_H2_PBE_G0W0_CHOLESKY_OFF.inp b/tests/QS/regtest-bse/TDA_H2_PBE_G0W0_CHOLESKY_OFF.inp new file mode 100644 index 0000000000..fe406f4abb --- /dev/null +++ b/tests/QS/regtest-bse/TDA_H2_PBE_G0W0_CHOLESKY_OFF.inp @@ -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 diff --git a/tests/QS/regtest-bse/TEST_FILES.toml b/tests/QS/regtest-bse/TEST_FILES.toml index e303207e44..095e4fcbc1 100644 --- a/tests/QS/regtest-bse/TEST_FILES.toml +++ b/tests/QS/regtest-bse/TEST_FILES.toml @@ -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