diff --git a/src/almo_scf_methods.F b/src/almo_scf_methods.F index fd80a75780..1aa765ad17 100644 --- a/src/almo_scf_methods.F +++ b/src/almo_scf_methods.F @@ -2216,9 +2216,8 @@ CONTAINS index1_start, n_domain_mos, naos, & nblkrows_tot, ndomains, neighbor, row INTEGER, DIMENSION(:), POINTER :: nmos - LOGICAL :: matrix_r_required, & - matrix_s_inv_required, & - matrix_trimmer_required, my_use_trimmer + LOGICAL :: eps_zero_eigenvalues_required, matrix_r_required, matrix_s_half_required, & + matrix_s_inv_half_required, matrix_s_inv_required, matrix_trimmer_required, my_use_trimmer REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Minv, proj_array TYPE(domain_submatrix_type), ALLOCATABLE, & DIMENSION(:) :: subm_main, subm_tmp, subm_tmp2 @@ -2230,6 +2229,9 @@ CONTAINS my_use_trimmer = use_trimmer ENDIF + matrix_s_inv_half_required = .FALSE. + matrix_s_half_required = .FALSE. + eps_zero_eigenvalues_required = .FALSE. matrix_s_inv_required = .FALSE. matrix_trimmer_required = .FALSE. matrix_r_required = .FALSE. @@ -2240,7 +2242,23 @@ CONTAINS matrix_trimmer_required = .TRUE. CPABORT("TRIMMED PRECONDITIONER DISABLED!") ENDIF + ! tie the following optional arguments together to prevent bad calls + IF (PRESENT(bad_modes_projector_down)) THEN + matrix_s_inv_half_required = .TRUE. + matrix_s_half_required = .TRUE. + eps_zero_eigenvalues_required = .TRUE. + ENDIF + ! check if all required optional arguments are provided + IF (.NOT. PRESENT(subm_s_inv_half) .AND. matrix_s_inv_half_required) THEN + CPABORT("S_inv_half SUBMATRICES ARE REQUIRED") + ENDIF + IF (.NOT. PRESENT(subm_s_half) .AND. matrix_s_half_required) THEN + CPABORT("S_half SUBMATRICES ARE REQUIRED") + ENDIF + IF (.NOT. PRESENT(eps_zero_eigenvalues) .AND. eps_zero_eigenvalues_required) THEN + CPABORT("EPS_ZERO_EIGENVALUES IS REQUIRED") + ENDIF IF (.NOT. PRESENT(subm_s_inv) .AND. matrix_s_inv_required) THEN CPABORT("S_inv SUBMATRICES ARE REQUIRED") ENDIF @@ -2347,50 +2365,26 @@ CONTAINS IF (PRESENT(bad_modes_projector_down)) THEN ALLOCATE (proj_array(naos, naos)) - ENDIF - - IF (PRESENT(eps_zero_eigenvalues)) THEN - IF (PRESENT(subm_s_inv_half)) THEN - CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, & - range1=nmos(idomain), range2=n_domain_mos, & - range1_thr=eps_zero_eigenvalues, & - bad_modes_projector_down=proj_array, & - s_inv_half=subm_s_inv_half(idomain)%mdata, & - s_half=subm_s_half(idomain)%mdata & - !metric_inv=subm_s_inv(idomain)%mdata & - ) - ELSE - CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, & - range1=nmos(idomain), range2=n_domain_mos, & - range1_thr=eps_zero_eigenvalues, & - bad_modes_projector_down=proj_array & - ) - ENDIF + CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, & + range1=nmos(idomain), range2=n_domain_mos, & + range1_thr=eps_zero_eigenvalues, & + bad_modes_projector_down=proj_array, & + s_inv_half=subm_s_inv_half(idomain)%mdata, & + s_half=subm_s_half(idomain)%mdata & + ) ELSE - IF (PRESENT(subm_s_inv_half)) THEN - CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, & - range1=nmos(idomain), range2=n_domain_mos, & - bad_modes_projector_down=proj_array, & - s_inv_half=subm_s_inv_half(idomain)%mdata, & - s_half=subm_s_inv(idomain)%mdata) - ELSE - CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, & - range1=nmos(idomain), range2=n_domain_mos, & - bad_modes_projector_down=proj_array) - ENDIF - ENDIF ! eps_zero_eigenvalues + CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, & + range1=nmos(idomain), range2=n_domain_mos) + ENDIF !!!TRIM ENDIF CALL copy_submatrices(subm_main(idomain), preconditioner(idomain), .FALSE.) CALL copy_submatrix_data(Minv, preconditioner(idomain)) + DEALLOCATE (Minv) IF (PRESENT(bad_modes_projector_down)) THEN CALL copy_submatrices(subm_main(idomain), bad_modes_projector_down(idomain), .FALSE.) CALL copy_submatrix_data(proj_array, bad_modes_projector_down(idomain)) - ENDIF - - DEALLOCATE (Minv) - IF (PRESENT(bad_modes_projector_down)) THEN DEALLOCATE (proj_array) ENDIF @@ -2409,8 +2403,6 @@ CONTAINS ! CALL dbcsr_release(matrix_r) !ENDIF - !RZK-warning do we need a barrier here ? - CALL timestop(handle) END SUBROUTINE construct_domain_preconditioner @@ -2760,8 +2752,8 @@ CONTAINS INTEGER, INTENT(IN) :: N, method INTEGER, INTENT(IN), OPTIONAL :: range1, range2 REAL(KIND=dp), INTENT(IN), OPTIONAL :: range1_thr, shift - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & - INTENT(INOUT), OPTIONAL :: bad_modes_projector_down + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & + OPTIONAL :: bad_modes_projector_down REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & OPTIONAL :: s_inv_half, s_half @@ -2770,7 +2762,7 @@ CONTAINS INTEGER :: handle, ii, INFO, jj, LWORK, range1_eiv, & range2_eiv, range3_eiv, unit_nr - LOGICAL :: use_both, use_ranges + LOGICAL :: use_both, use_ranges_only, use_thr_only REAL(KIND=dp) :: my_shift REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, WORK REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: temp1, temp2, temp3, temp4 @@ -2787,15 +2779,30 @@ CONTAINS END IF IF (method .EQ. 1) THEN + + IF ((PRESENT(range1) .AND. (.NOT. PRESENT(range2))) .OR. (PRESENT(range2) .AND. (.NOT. PRESENT(range1)))) THEN + CPABORT("range1 and range2 must be provided together") + ENDIF + IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN use_both = .TRUE. + use_thr_only = .FALSE. + use_ranges_only = .FALSE. ELSE use_both = .FALSE. + IF (PRESENT(range1)) THEN - use_ranges = .TRUE. + use_ranges_only = .TRUE. ELSE - use_ranges = .FALSE. + use_ranges_only = .FALSE. ENDIF + + IF (PRESENT(range1_thr)) THEN + use_thr_only = .TRUE. + ELSE + use_thr_only = .FALSE. + ENDIF + ENDIF IF ((PRESENT(s_half) .AND. (.NOT. PRESENT(s_inv_half))) .OR. (PRESENT(s_inv_half) .AND. (.NOT. PRESENT(s_half)))) THEN @@ -2812,7 +2819,7 @@ CONTAINS INFO = 0 SELECT CASE (method) - CASE (0) + CASE (0) ! Inversion via cholesky factorization CALL DPOTRF('L', N, Ainv, N, INFO) IF (INFO .NE. 0) THEN @@ -2863,9 +2870,7 @@ CONTAINS ! invert eigenvalues and use eigenvectors to compute pseudo Ainv ! project out near-zero eigenvalue modes ALLOCATE (temp2(N, N)) - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) ALLOCATE (temp3(N, N)) - ENDIF + IF (PRESENT(bad_modes_projector_down)) ALLOCATE (temp3(N, N)) temp2(1:N, 1:N) = Ainv(1:N, 1:N) range1_eiv = 0 @@ -2876,70 +2881,55 @@ CONTAINS DO jj = 1, N IF ((jj .LE. range2) .AND. (eigenvalues(jj) .LT. range1_thr)) THEN temp1(jj, :) = temp2(:, jj)*0.0_dp - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp - ENDIF + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp range1_eiv = range1_eiv+1 ELSE temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj)+my_shift) - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp - ENDIF + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp range2_eiv = range2_eiv+1 ENDIF ENDDO ELSE - IF (use_ranges) THEN + IF (use_ranges_only) THEN DO jj = 1, N IF (jj .LE. range1) THEN temp1(jj, :) = temp2(:, jj)*0.0_dp - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp - ENDIF + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp range1_eiv = range1_eiv+1 ELSE IF (jj .LE. range2) THEN temp1(jj, :) = temp2(:, jj)*1.0_dp - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp - ENDIF + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp range2_eiv = range2_eiv+1 ELSE temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj)+my_shift) - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp - ENDIF + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp range3_eiv = range3_eiv+1 ENDIF ENDDO - ELSE + ELSE IF (use_thr_only) THEN DO jj = 1, N IF (eigenvalues(jj) .LT. range1_thr) THEN temp1(jj, :) = temp2(:, jj)*0.0_dp - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp - ENDIF + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp range1_eiv = range1_eiv+1 ELSE temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj)+my_shift) - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp - ENDIF - + IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp range2_eiv = range2_eiv+1 ENDIF ENDDO + ELSE ! no ranges, no thresholds + CPABORT("Invert using Cholesky. It would be faster.") ENDIF ENDIF !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) THEN - IF (PRESENT(s_half)) THEN - CALL DSYMM('L', 'U', N, N, 1.0_dp, s_half, N, temp2, N, 0.0_dp, Ainv, N) - CALL DSYMM('R', 'U', N, N, 1.0_dp, s_half, N, temp3, N, 0.0_dp, temp4, N) - CALL DGEMM('N', 'N', N, N, N, 1.0_dp, Ainv, N, temp4, N, 0.0_dp, bad_modes_projector_down, N) - ELSE - CALL DGEMM('N', 'N', N, N, N, 1.0_dp, temp2, N, temp3, N, 0.0_dp, bad_modes_projector_down, N) - ENDIF + IF (PRESENT(s_half)) THEN + CALL DSYMM('L', 'U', N, N, 1.0_dp, s_half, N, temp2, N, 0.0_dp, Ainv, N) + CALL DSYMM('R', 'U', N, N, 1.0_dp, s_half, N, temp3, N, 0.0_dp, temp4, N) + CALL DGEMM('N', 'N', N, N, N, 1.0_dp, Ainv, N, temp4, N, 0.0_dp, bad_modes_projector_down, N) + ELSE + CALL DGEMM('N', 'N', N, N, N, 1.0_dp, temp2, N, temp3, N, 0.0_dp, bad_modes_projector_down, N) ENDIF ENDIF @@ -2951,9 +2941,7 @@ CONTAINS CALL DGEMM('N', 'N', N, N, N, 1.0_dp, temp2, N, temp1, N, 0.0_dp, Ainv, N) ENDIF DEALLOCATE (temp1, temp2, temp4) - IF (PRESENT(bad_modes_projector_down)) THEN - IF (ALLOCATED(bad_modes_projector_down)) DEALLOCATE (temp3) - ENDIF + IF (PRESENT(bad_modes_projector_down)) DEALLOCATE (temp3) DEALLOCATE (eigenvalues) CASE DEFAULT diff --git a/src/almo_scf_optimizer.F b/src/almo_scf_optimizer.F index 65ab3feed5..52f2072b62 100644 --- a/src/almo_scf_optimizer.F +++ b/src/almo_scf_optimizer.F @@ -906,6 +906,13 @@ CONTAINS nspins = almo_scf_env%nspins + ! if unprojected XALMOs are optimized + ! then we must use the "blissful_neglect" procedure + blissful_neglect = .FALSE. + IF (my_special_case .EQ. xalmo_case_normal .AND. .NOT. assume_t0_q0x) THEN + blissful_neglect = .TRUE. + ENDIF + IF (unit_nr > 0) THEN WRITE (unit_nr, *) SELECT CASE (my_special_case) @@ -916,8 +923,13 @@ CONTAINS WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 20), & " Optimization of fully delocalized MOs ", REPEAT("-", 20) CASE (xalmo_case_normal) - WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 27), & - " Optimization of XALMOs ", REPEAT("-", 28) + IF (blissful_neglect) THEN + WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 25), & + " LCP optimization of XALMOs ", REPEAT("-", 26) + ELSE + WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 27), & + " Optimization of XALMOs ", REPEAT("-", 28) + ENDIF END SELECT WRITE (unit_nr, *) WRITE (unit_nr, '(T2,A13,A6,A23,A14,A14,A9)') "Method", "Iter", & @@ -930,13 +942,6 @@ CONTAINS optimize_theta = almo_scf_env%logical05 eps_skip_gradients = almo_scf_env%real01 - ! if unprojected XALMOs are optimized then compute both - ! then we must use the "blissful_neglect" procedure - blissful_neglect = .FALSE. - IF (my_special_case .EQ. xalmo_case_normal .AND. .NOT. assume_t0_q0x) THEN - blissful_neglect = .TRUE. - ENDIF - ! penalty amplitude adjusts the strenght of volume conservation ! the following guidelines are useful ! A = T for n = 2