From e19d499b6491daad1212c568ee67f1490bf079c2 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Thu, 9 Mar 2023 13:57:14 +0100 Subject: [PATCH] Refactoring of linres solver (#2665) --- src/mp2_cphf.F | 8 ++- src/preconditioner.F | 84 +++++++++++++++++++------- src/preconditioner_makes.F | 6 +- src/qs_density_matrices.F | 35 ++++++++--- src/qs_linres_methods.F | 68 ++++++++++++++------- src/qs_scf_diagonalization.F | 3 +- src/response_solver.F | 49 ++++++++++----- src/rpa_im_time_force_methods.F | 5 +- tests/QS/regtest-mp2-grad-2/TEST_FILES | 2 +- tests/TEST_DIRS | 10 +-- 10 files changed, 187 insertions(+), 83 deletions(-) diff --git a/src/mp2_cphf.F b/src/mp2_cphf.F index 712ed7fe68..bef5fd6b14 100644 --- a/src/mp2_cphf.F +++ b/src/mp2_cphf.F @@ -112,7 +112,8 @@ MODULE mp2_cphf USE qs_ks_types, ONLY: qs_ks_env_type,& set_ks_env USE qs_linres_types, ONLY: linres_control_type - USE qs_mo_types, ONLY: mo_set_type + USE qs_mo_types, ONLY: get_mo_set,& + mo_set_type USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type USE qs_overlap, ONLY: build_overlap_matrix USE qs_p_env_methods, ONLY: p_env_check_i_alloc,& @@ -1316,7 +1317,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'update_mp2_forces' INTEGER :: alpha, beta, handle, idir, iounit, & - ispin, nspins + ispin, nocc, nspins INTEGER, DIMENSION(3) :: comp LOGICAL :: do_exx, do_hfx, use_virial REAL(KIND=dp) :: e_dummy, e_hartree, e_xc, ehartree, exc @@ -1948,8 +1949,9 @@ CONTAINS ! Finish matrix_w_mp2 with occ-occ block DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), homo=nocc, nmo=alpha) CALL calculate_whz_matrix(mos(ispin)%mo_coeff, p_env%kpp1(ispin)%matrix, & - p_env%w1(1)%matrix, 1.0_dp) + p_env%w1(1)%matrix, 1.0_dp, nocc) END DO IF (debug_forces .AND. use_virial) e_dummy = third_tr(virial%pv_virial) diff --git a/src/preconditioner.F b/src/preconditioner.F index 6847d6de13..46fc2451aa 100644 --- a/src/preconditioner.F +++ b/src/preconditioner.F @@ -16,7 +16,13 @@ MODULE preconditioner USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm - USE cp_fm_types, ONLY: cp_fm_get_info,& + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type + USE cp_fm_types, ONLY: cp_fm_create,& + cp_fm_get_info,& + cp_fm_release,& + cp_fm_to_fm,& cp_fm_type USE dbcsr_api, ONLY: dbcsr_p_type,& dbcsr_type @@ -45,7 +51,8 @@ MODULE preconditioner qs_environment_type USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues USE qs_mo_types, ONLY: get_mo_set,& - mo_set_type + mo_set_type,& + set_mo_set #include "./base/base_uses.f90" IMPLICIT NONE @@ -106,20 +113,22 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner' - INTEGER :: handle, k, my_solver_type + INTEGER :: handle, k, my_solver_type, nao, nhomo LOGICAL :: my_convert_precond_to_dbcsr, & needs_full_spectrum, needs_homo, & use_mo_coeff_b REAL(KIND=dp) :: energy_homo - REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_ot + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot + TYPE(cp_fm_struct_type), POINTER :: fm_struct + TYPE(cp_fm_type) :: mo_occ TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(dbcsr_type), POINTER :: mo_coeff_b CALL timeset(routineN, handle) - CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b) + CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo) use_mo_coeff_b = mo_set%use_mo_coeff_b - CALL cp_fm_get_info(mo_coeff, ncol_global=k) + CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao) ! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix my_convert_precond_to_dbcsr = .FALSE. @@ -152,9 +161,9 @@ CONTAINS CPABORT("The preconditioner is unknown ...") END SELECT + ALLOCATE (eigenvalues_ot(k)) energy_homo = 0.0_dp IF (needs_full_spectrum) THEN - ALLOCATE (eigenvalues_ot(k)) ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either ! the case of non-equivalent MOs but also for the derivate ! we could already have all eigenvalues e.g. full_all and we could skip this @@ -170,7 +179,10 @@ CONTAINS CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, & eigenvalues_ot, do_rotation=.FALSE.) END IF - IF (k > 0) energy_homo = eigenvalues_ot(k) + IF (k > 0) THEN + CPASSERT(nhomo > 0 .AND. nhomo <= k) + energy_homo = eigenvalues_ot(nhomo) + END IF ELSE IF (needs_homo) THEN CPABORT("Not yet implemented") @@ -183,16 +195,24 @@ CONTAINS my_solver_type = solver_type preconditioner_env%in_use = precon_type preconditioner_env%cholesky_use = cholesky_reduce -!dbg -! write(*,*) PRESENT(chol_type) -! IF(PRESENT(chol_type)) write(*,*) chol_type, cholesky_reduce -! stop -!dbg IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type preconditioner_env%in_use = precon_type - CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, & - energy_homo, eigenvalues_ot, energy_gap, & - my_solver_type) + IF (nhomo == k) THEN + CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, & + energy_homo, eigenvalues_ot, energy_gap, my_solver_type) + ELSE + CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, & + context=preconditioner_env%ctxt, & + para_env=preconditioner_env%para_env) + CALL cp_fm_create(mo_occ, fm_struct) + CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo) + CALL cp_fm_struct_release(fm_struct) + ! + CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, & + energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type) + ! + CALL cp_fm_release(mo_occ) + END IF CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h) @@ -206,9 +226,7 @@ CONTAINS preconditioner_env%para_env, preconditioner_env%ctxt) END IF - IF (needs_full_spectrum) THEN - DEALLOCATE (eigenvalues_ot) - END IF + DEALLOCATE (eigenvalues_ot) CALL timestop(handle) @@ -282,14 +300,15 @@ CONTAINS !> \param has_unit_metric ... !> \param convert_to_dbcsr ... !> \param chol_type ... +!> \param full_mo_set ... ! ************************************************************************************************** SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, & ot_preconditioner, prec_type, solver_type, & energy_gap, nspins, has_unit_metric, & - convert_to_dbcsr, chol_type) + convert_to_dbcsr, chol_type, full_mo_set) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos + TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner INTEGER, INTENT(IN) :: prec_type, solver_type @@ -297,13 +316,16 @@ CONTAINS INTEGER, INTENT(IN) :: nspins LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr INTEGER, INTENT(IN), OPTIONAL :: chol_type + LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner' CHARACTER(LEN=default_string_length) :: msg INTEGER :: handle, icall, ispin, n_loops + INTEGER, DIMENSION(5) :: nocc, norb LOGICAL :: do_co_rotate, my_convert_to_dbcsr, & - my_has_unit_metric, use_mo_coeff_b + my_full_mo_set, my_has_unit_metric, & + use_mo_coeff_b TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic @@ -317,6 +339,8 @@ CONTAINS IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric my_convert_to_dbcsr = .TRUE. IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr + my_full_mo_set = .FALSE. + IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set CALL get_qs_env(qs_env, & dft_control=dft_control, & @@ -336,6 +360,15 @@ CONTAINS matrix_t => kinetic(1)%matrix END IF + ! use full set of MOs or just occupied MOs + nocc = 0 + norb = 0 + IF (my_full_mo_set) THEN + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin)) + CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin)) + END DO + END IF !determines how often make preconditioner is called, spin dependent methods have to be called twice n_loops = 1 IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins @@ -399,6 +432,13 @@ CONTAINS END DO END SELECT + ! reset homo values + IF (my_full_mo_set) THEN + DO ispin = 1, nspins + CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin)) + END DO + END IF + CALL timestop(handle) END SUBROUTINE prepare_preconditioner diff --git a/src/preconditioner_makes.F b/src/preconditioner_makes.F index 7274ab249c..7efb65f8ea 100644 --- a/src/preconditioner_makes.F +++ b/src/preconditioner_makes.F @@ -83,7 +83,7 @@ CONTAINS TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t TYPE(cp_fm_type), INTENT(IN) :: mo_coeff REAL(KIND=dp) :: energy_homo - REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_ot + REAL(KIND=dp), DIMENSION(:) :: eigenvalues_ot REAL(KIND=dp) :: energy_gap INTEGER :: my_solver_type @@ -377,7 +377,7 @@ CONTAINS TYPE(preconditioner_type) :: preconditioner_env TYPE(cp_fm_type), INTENT(IN) :: matrix_c0 TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s - REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals + REAL(KIND=dp), DIMENSION(:) :: c0_evals REAL(KIND=dp) :: energy_gap CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all' @@ -552,7 +552,7 @@ CONTAINS TYPE(preconditioner_type) :: preconditioner_env TYPE(cp_fm_type), INTENT(IN) :: matrix_c0 TYPE(dbcsr_type), POINTER :: matrix_h - REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals + REAL(KIND=dp), DIMENSION(:) :: c0_evals REAL(KIND=dp) :: energy_gap CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho' diff --git a/src/qs_density_matrices.F b/src/qs_density_matrices.F index 7437029b78..7bf8f3325c 100644 --- a/src/qs_density_matrices.F +++ b/src/qs_density_matrices.F @@ -368,25 +368,40 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'calculate_wz_matrix' - INTEGER :: handle, ncol, nrow + INTEGER :: handle, ncol, nocc, nrow TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp TYPE(cp_fm_type) :: ksmat, scrv CALL timeset(routineN, handle) +! CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow) +! CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors") +! CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, & +! para_env=mo_set%mo_coeff%matrix_struct%para_env, & +! context=mo_set%mo_coeff%matrix_struct%context) +! CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS") +! CALL cp_fm_struct_release(fm_struct_tmp) +! CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol) +! CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat) +! CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv) +! CALL dbcsr_set(w_matrix, 0.0_dp) +! CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, & +! ncol=mo_set%homo, symmetry_mode=1) +! CALL cp_fm_release(scrv) +! CALL cp_fm_release(ksmat) CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow) + nocc = mo_set%homo CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors") - CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, & + CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nocc, ncol_global=nocc, & para_env=mo_set%mo_coeff%matrix_struct%para_env, & context=mo_set%mo_coeff%matrix_struct%context) CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS") CALL cp_fm_struct_release(fm_struct_tmp) - CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol) - CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat) - CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv) + CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc) + CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat) + CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv) CALL dbcsr_set(w_matrix, 0.0_dp) - CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, & - ncol=mo_set%homo, symmetry_mode=1) + CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1) CALL cp_fm_release(scrv) CALL cp_fm_release(ksmat) @@ -401,15 +416,17 @@ CONTAINS !> \param hzm ... !> \param w_matrix sparse matrix !> \param focc ... +!> \param nocc ... !> \par History !> adapted from calculate_w_matrix_1 !> \author JGH ! ************************************************************************************************** - SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc) + SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc) TYPE(cp_fm_type), INTENT(IN) :: c0vec TYPE(dbcsr_type), POINTER :: hzm, w_matrix REAL(KIND=dp), INTENT(IN) :: focc + INTEGER, INTENT(IN) :: nocc CHARACTER(len=*), PARAMETER :: routineN = 'calculate_whz_matrix' @@ -424,6 +441,8 @@ CONTAINS CALL cp_fm_create(hcvec, c0vec%matrix_struct, "hcvec") CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb) + CPASSERT(nocc <= norb .AND. nocc > 0) + norb = nocc CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, & ncol_global=norb, para_env=fm_struct%para_env) CALL cp_fm_create(chcmat, fm_struct_mat) diff --git a/src/qs_linres_methods.F b/src/qs_linres_methods.F index 8c2029a237..26cb2b8feb 100644 --- a/src/qs_linres_methods.F +++ b/src/qs_linres_methods.F @@ -199,6 +199,9 @@ CONTAINS !> \brief scf loop to optimize the first order wavefunctions (psi1) !> given a perturbation as an operator applied to the ground !> state orbitals (h1_psi0) +!> psi1 is defined wrt psi0_order (can be a subset of the occupied space) +!> The reference ground state is defined through qs_env (density and ground state MOs) +!> psi1 is orthogonal to the occupied orbitals in the ground state MO set (qs_env%mos) !> \param p_env ... !> \param qs_env ... !> \param psi1 ... @@ -219,8 +222,8 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_solver' - INTEGER :: handle, ispin, iter, maxnmo, maxnmo_o, & - nao, ncol, nmo, nspins + INTEGER :: handle, ispin, iter, maxnmo, nao, ncol, & + nmo, nocc, nspins LOGICAL :: restart REAL(dp) :: alpha, beta, norm_res, t1, t2 REAL(dp), DIMENSION(:), POINTER :: tr_pAp, tr_rz0, tr_rz00, tr_rz1 @@ -253,25 +256,18 @@ CONTAINS mos=mos) nspins = dft_control%nspins - CALL get_mo_set(mos(1), nao=nao) + CALL cp_fm_get_info(psi1(1), nrow_global=nao) maxnmo = 0 - maxnmo_o = 0 DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), nmo=ncol) + CALL cp_fm_get_info(psi1(ispin), ncol_global=ncol) maxnmo = MAX(maxnmo, ncol) - CALL cp_fm_get_info(psi0_order(ispin), ncol_global=ncol) - maxnmo_o = MAX(maxnmo_o, ncol) END DO ! CALL check_p_env_init(p_env, linres_control, nspins) ! ! allocate the vectors ALLOCATE (tr_pAp(nspins), tr_rz0(nspins), tr_rz00(nspins), tr_rz1(nspins), & - r(nspins), p(nspins), z(nspins), Ap(nspins), mo_coeff_array(nspins)) - DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff) - mo_coeff_array(ispin) = mo_coeff - END DO + r(nspins), p(nspins), z(nspins), Ap(nspins)) ! DO ispin = 1, nspins CALL cp_fm_create(r(ispin), psi1(ispin)%matrix_struct) @@ -280,11 +276,24 @@ CONTAINS CALL cp_fm_create(Ap(ispin), psi1(ispin)%matrix_struct) END DO ! - ! compute S*C0, C0_order'*H*C0_order (this should be done once for all) - ALLOCATE (chc(nspins), Sc(nspins)) + ! build C0 occupied vectors and S*C0 matrix + ALLOCATE (Sc(nspins), mo_coeff_array(nspins)) DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo) - CALL cp_fm_create(Sc(ispin), mo_coeff%matrix_struct) + CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc) + NULLIFY (tmp_fm_struct) + CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, & + ncol_global=nocc, para_env=para_env, & + context=mo_coeff%matrix_struct%context) + CALL cp_fm_create(mo_coeff_array(ispin), tmp_fm_struct) + CALL cp_fm_to_fm(mo_coeff, mo_coeff_array(ispin), nocc) + CALL cp_fm_create(Sc(ispin), tmp_fm_struct) + CALL cp_fm_struct_release(tmp_fm_struct) + END DO + ! + ! Allocate C0_order'*H*C0_order + ALLOCATE (chc(nspins)) + DO ispin = 1, nspins + CALL cp_fm_get_info(psi1(ispin), ncol_global=nmo) NULLIFY (tmp_fm_struct) CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, & ncol_global=nmo, para_env=para_env, & @@ -305,9 +314,8 @@ CONTAINS END ASSOCIATE ! ! S * C0 - CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff) - CALL cp_fm_get_info(mo_coeff, ncol_global=ncol) - CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, Sc(ispin), ncol) + CALL cp_fm_get_info(mo_coeff_array(ispin), ncol_global=ncol) + CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff_array(ispin), Sc(ispin), ncol) END DO ! ! header @@ -384,7 +392,7 @@ CONTAINS CALL cp_fm_trace(r(ispin), z(ispin), tr_rz0(ispin)) END DO IF (SUM(tr_rz0) < 0.0_dp) CPABORT("tr(r_j*z_j) < 0") - norm_res = ABS(SUM(tr_rz0))/SQRT(REAL(nspins*nao*maxnmo_o, dp)) + norm_res = ABS(SUM(tr_rz0))/SQRT(REAL(nspins*nao*maxnmo, dp)) ! alpha = 0.0_dp restart = .FALSE. @@ -463,7 +471,6 @@ CONTAINS CALL apply_op(qs_env, p_env, psi0_order, psi1, Ap, chc) ! DO ispin = 1, nspins - CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff) CALL cp_fm_to_fm(h1_psi0(ispin), r(ispin)) CALL cp_fm_scale_and_add(-1.0_dp, r(ispin), -1.0_dp, Ap(ispin)) END DO @@ -506,7 +513,7 @@ CONTAINS CALL cp_fm_trace(r(ispin), z(ispin), tr_rz1(ispin)) END DO IF (SUM(tr_rz1) < 0.0_dp) CPABORT("tr(r_j+1*z_j+1) < 0") - norm_res = SUM(tr_rz1)/SQRT(REAL(nspins*nao*maxnmo_o, dp)) + norm_res = SUM(tr_rz1)/SQRT(REAL(nspins*nao*maxnmo, dp)) ! ! beta = tr(r_j+1*z_j+1) / tr(r_j*z_j) IF (SUM(tr_rz0) < 1.0e-10_dp) THEN @@ -539,6 +546,7 @@ CONTAINS CALL cp_fm_release(z(ispin)) CALL cp_fm_release(Ap(ispin)) ! + CALL cp_fm_release(mo_coeff_array(ispin)) CALL cp_fm_release(Sc(ispin)) CALL cp_fm_release(chc(ispin)) END DO @@ -565,7 +573,8 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op' - INTEGER :: handle, ispin, nspins + INTEGER :: handle, ispin, nc1, nc2, nc3, nc4, nr1, & + nr2, nr3, nr4, nspins REAL(dp) :: chksum TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(dft_control_type), POINTER :: dft_control @@ -583,6 +592,19 @@ CONTAINS nspins = dft_control%nspins + DO ispin = 1, nspins + !c0, v, Av, chc + CALL cp_fm_get_info(c0(ispin), ncol_global=nc1, nrow_global=nr1) + CALL cp_fm_get_info(v(ispin), ncol_global=nc2, nrow_global=nr2) + CALL cp_fm_get_info(Av(ispin), ncol_global=nc3, nrow_global=nr3) + CALL cp_fm_get_info(chc(ispin), ncol_global=nc4, nrow_global=nr4) + IF (.NOT. (nc1 == nc2 .AND. nr1 == nr2 .AND. nc1 == nc3 .AND. nr1 == nr3 & + .AND. nc4 == nr4 .AND. nc1 <= nc4)) THEN + CALL cp_abort(__LOCATION__, & + "Number of vectors inconsistent or CHC matrix too small") + END IF + END DO + ! apply the uncoupled operator DO ispin = 1, nspins CALL apply_op_1(v(ispin), Av(ispin), matrix_ks(ispin)%matrix, & diff --git a/src/qs_scf_diagonalization.F b/src/qs_scf_diagonalization.F index 2e9cc30114..2839343d6c 100644 --- a/src/qs_scf_diagonalization.F +++ b/src/qs_scf_diagonalization.F @@ -1641,7 +1641,8 @@ CONTAINS scf_env%block_davidson_env(1)%prec_type, & scf_env%block_davidson_env(1)%solver_type, & scf_env%block_davidson_env(1)%energy_gap, nspins, & - convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos) + convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, & + full_mo_set=.TRUE.) END IF DO ispin = 1, nspins diff --git a/src/response_solver.F b/src/response_solver.F index 2e18bb6c9d..25d3cf4e26 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -32,12 +32,15 @@ MODULE response_solver cp_dbcsr_sm_fm_multiply,& dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set - USE cp_fm_struct, ONLY: cp_fm_struct_type + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type USE cp_fm_types, ONLY: cp_fm_create,& cp_fm_get_info,& cp_fm_init_random,& cp_fm_release,& cp_fm_set_all,& + cp_fm_to_fm,& cp_fm_type USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_unit_nr,& @@ -546,7 +549,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'response_equation_new' - INTEGER :: handle, ispin, nao, nao_aux, nspins + INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins LOGICAL :: should_stop TYPE(admm_type), POINTER :: admm_env TYPE(cp_fm_struct_type), POINTER :: fm_struct @@ -555,13 +558,14 @@ CONTAINS TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(dft_control_type), POINTER :: dft_control TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(mp_para_env_type), POINTER :: para_env CALL timeset(routineN, handle) NULLIFY (dft_control, matrix_ks, mo_coeff, mos) CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks, & - matrix_s=matrix_s, mos=mos) + para_env=para_env, matrix_s=matrix_s, mos=mos) nspins = dft_control%nspins ! Initialize vectors: @@ -569,11 +573,16 @@ CONTAINS ! psi1 : The "perturbed" linear response orbitals ALLOCATE (psi0(nspins), psi1(nspins)) DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff) - psi0(ispin) = mo_coeff - CALL cp_fm_get_info(mo_coeff, matrix_struct=fm_struct) + CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc) + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, nrow_global=nao, & + ncol_global=nocc, para_env=para_env, & + context=mo_coeff%matrix_struct%context) + CALL cp_fm_create(psi0(ispin), fm_struct) + CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc) CALL cp_fm_create(psi1(ispin), fm_struct) CALL cp_fm_set_all(psi1(ispin), 0.0_dp) + CALL cp_fm_struct_release(fm_struct) END DO should_stop = .FALSE. @@ -617,6 +626,7 @@ CONTAINS DO ispin = 1, nspins CALL cp_fm_release(cpmos(ispin)) CALL cp_fm_release(psi1(ispin)) + CALL cp_fm_release(psi0(ispin)) END DO DEALLOCATE (psi0, psi1) @@ -643,7 +653,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'response_equation' - INTEGER :: handle, ispin, nao, nao_aux, nspins + INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins LOGICAL :: should_stop TYPE(admm_type), POINTER :: admm_env TYPE(cp_fm_struct_type), POINTER :: fm_struct @@ -653,6 +663,7 @@ CONTAINS TYPE(dft_control_type), POINTER :: dft_control TYPE(linres_control_type), POINTER :: linres_control TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb @@ -688,7 +699,7 @@ CONTAINS CALL p_env_psi0_changed(p_env, qs_env) p_env%new_preconditioner = .TRUE. - CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos) + CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, mos=mos) ! nspins = dft_control%nspins @@ -697,11 +708,16 @@ CONTAINS ! psi1 : The "perturbed" linear response orbitals ALLOCATE (psi0(nspins), psi1(nspins)) DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff) - psi0(ispin) = mo_coeff - CALL cp_fm_get_info(mo_coeff, matrix_struct=fm_struct) + CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc) + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, nrow_global=nao, & + ncol_global=nocc, para_env=para_env, & + context=mo_coeff%matrix_struct%context) + CALL cp_fm_create(psi0(ispin), fm_struct) + CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc) CALL cp_fm_create(psi1(ispin), fm_struct) CALL cp_fm_set_all(psi1(ispin), 0.0_dp) + CALL cp_fm_struct_release(fm_struct) END DO should_stop = .FALSE. @@ -765,6 +781,7 @@ CONTAINS p_env%w1(ispin)%matrix) END DO DO ispin = 1, nspins + CALL cp_fm_release(psi0(ispin)) CALL cp_fm_release(psi1(ispin)) END DO DEALLOCATE (psi0, psi1) @@ -811,7 +828,7 @@ CONTAINS INTEGER :: handle, iounit, ispin, mspin, n_rep_hf, & nao, nao_aux, natom, nder, nimages, & - nspins + nocc, nspins INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: debug_forces, debug_stress, distribute_fock_matrix, do_ex, do_hfx, gapw, gapw_xc, & hfx_treat_lsd_in_core, resp_only, s_mstruct_changed, use_virial @@ -2036,8 +2053,9 @@ CONTAINS IF (nspins == 1) focc = 2.0_dp CALL get_qs_env(qs_env, mos=mos) DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), homo=nocc) CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, & - matrix_wz(ispin)%matrix, focc) + matrix_wz(ispin)%matrix, focc, nocc) END DO END IF IF (nspins == 2) THEN @@ -2110,7 +2128,7 @@ CONTAINS INTEGER :: atom_a, handle, iatom, ikind, iounit, & is, ispin, na, natom, natorb, nimages, & - nkind, ns, nsgf, nspins + nkind, nocc, ns, nsgf, nspins INTEGER, DIMENSION(25) :: lao INTEGER, DIMENSION(5) :: occ LOGICAL :: debug_forces, do_ex, use_virial @@ -2259,8 +2277,9 @@ CONTAINS IF (nspins == 1) focc = 1.0_dp CALL get_qs_env(qs_env, mos=mos) DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), homo=nocc) CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, & - matrix_wz(ispin)%matrix, focc) + matrix_wz(ispin)%matrix, focc, nocc) END DO IF (nspins == 2) THEN CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, & diff --git a/src/rpa_im_time_force_methods.F b/src/rpa_im_time_force_methods.F index 9ee122b29d..7578b0bcb5 100644 --- a/src/rpa_im_time_force_methods.F +++ b/src/rpa_im_time_force_methods.F @@ -3003,7 +3003,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'update_im_time_forces' INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, & - nao_aux, nder, nimages, nspins + nao_aux, nder, nimages, nocc, nspins INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, & use_virial @@ -3639,8 +3639,9 @@ CONTAINS IF (nspins == 2) focc = 1.0_dp CALL get_qs_env(qs_env, mos=mos) DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), homo=nocc) CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, & - p_env%w1(ispin)%matrix, focc) + p_env%w1(ispin)%matrix, focc, nocc) END DO IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, 1.0_dp) diff --git a/tests/QS/regtest-mp2-grad-2/TEST_FILES b/tests/QS/regtest-mp2-grad-2/TEST_FILES index 149db6a623..092ae73662 100644 --- a/tests/QS/regtest-mp2-grad-2/TEST_FILES +++ b/tests/QS/regtest-mp2-grad-2/TEST_FILES @@ -2,7 +2,7 @@ O2_dyn.inp 11 2e-11 - O2_dyn_mme.inp 72 1e-8 0.50585576 CH_dyn_screen.inp 11 2e-04 -5.700371305153854 CH_dyn_screen_single_group.inp 11 2e-04 -5.700371305153854 -MOM_MP2_geoopt.inp 11 4e-06 -13.969926262647746 +MOM_MP2_geoopt.inp 11 4e-06 -13.969416309750072 H2O_MD_mme.inp 11 1e-10 -17.056807598425291 O2_dyn_ri-hfx.inp 11 1e-10 -31.518511363027905 #EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 6d37680a27..edfa09d3a2 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -120,6 +120,8 @@ QS/regtest-gpw-3 QS/regtest-ps-implicit-1-3 fftw3 QS/regtest-ps-implicit-2-2 QS/regtest-rtp-3 +QS/regtest-rtp-4 +QS/regtest-rtp-5 QS/regtest-pao-1 QMMM/QS/regtest-1 QS/regtest-ps-implicit-2-1 @@ -240,6 +242,8 @@ QS/regtest-hybrid-2 libint QS/regtest-chi-1 QS/regtest-cdft-4-2 DFTB/regtest-nonscc +QS/regtest-as libint +QS/regtest-as-dft libint QS/regtest-double-hybrid-3 libint QS/regtest-dft-filtermat SE/regtest-3-3 @@ -285,8 +289,6 @@ QS/regtest-hfx-2 libint QS/regtest-nmr-3 QS/regtest-hybrid-1 libint QS/regtest-rtp-1 -QS/regtest-as libint -QS/regtest-as-dft libint QS/regtest-tddfpt-prop QS/regtest-sym-3 QS/regtest-gapw-ext @@ -297,6 +299,7 @@ QS/regtest-sym-1 QS/regtest-sym-2 QS/regtest-ot-refine-2 NEB/regtest-2 +QS/regtest-lsroks libint QS/regtest-nmr-2 QS/regtest-cdft-3-1 QS/regtest-optbas @@ -320,7 +323,6 @@ QS/regtest-pao-3 DFTB/regtest-debug optimize_input/regtest-1 FE/regtest-2 -QS/regtest-lsroks libint NEB/regtest-3 QS/regtest-linearscaling SE/regtest-3-1 @@ -335,7 +337,5 @@ Fist/regtest-plumed2 plumed2 Fist/regtest-quip quip Fist/regtest-nequip libtorch Fist/regtest-16 -QS/regtest-rtp-4 -QS/regtest-rtp-5 DFTB/regtest-vdw QS/regtest-tddfpt-force-gapw