Refactoring of linres solver (#2665)

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 @@ SUBROUTINE update_mp2_forces(qs_env)
       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 @@ SUBROUTINE update_mp2_forces(qs_env)
 
       ! 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)

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 @@ SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, mat
 
       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 @@ SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, mat
          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 @@ SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, mat
             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 @@ SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, mat
       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 @@ SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, mat
                                    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,28 +300,32 @@ END SUBROUTINE restart_preconditioner
 !> \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
       REAL(dp), INTENT(IN)                               :: energy_gap
       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 @@ SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
       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 @@ SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
          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 @@ SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
          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

src/preconditioner_makes.F

@@ -83,7 +83,7 @@ SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, ma
       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 @@ SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_e
       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 @@ SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals
       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'

src/qs_density_matrices.F

@@ -368,25 +368,40 @@ SUBROUTINE calculate_wz_matrix(mo_set, psi1, ks_matrix, w_matrix)
 
       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 @@ END SUBROUTINE calculate_wz_matrix
 !> \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 @@ SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc)
 
       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)