diff --git a/src/almo_scf.F b/src/almo_scf.F index 3cf26c2531..a1d6906efa 100644 --- a/src/almo_scf.F +++ b/src/almo_scf.F @@ -5,7 +5,7 @@ ! ************************************************************************************************** !> \brief Routines for all ALMO-based SCF methods -!> RZK-warning marks unresolved issues +!> 'RZK-warning' marks unresolved issues !> \par History !> 2011.05 created [Rustam Z Khaliullin] !> \author Rustam Z Khaliullin @@ -17,6 +17,7 @@ MODULE almo_scf distribute_domains,& orthogonalize_mos USE almo_scf_optimizer, ONLY: almo_scf_block_diagonal,& + almo_scf_construct_nlmos,& almo_scf_xalmo_eigensolver,& almo_scf_xalmo_pcg,& almo_scf_xalmo_trustr @@ -44,6 +45,7 @@ MODULE almo_scf USE cp_blacs_env, ONLY: cp_blacs_env_release,& cp_blacs_env_retain USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm USE cp_fm_types, ONLY: cp_fm_type USE cp_log_handling, ONLY: cp_get_default_logger,& @@ -53,9 +55,11 @@ MODULE almo_scf cp_para_env_retain USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: & - dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, & + dbcsr_add, dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, & dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, & - dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_multiply, dbcsr_nblkcols_total, & + dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_init_random, & + dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, & + dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_nblkcols_total, & dbcsr_nblkrows_total, dbcsr_p_type, dbcsr_release, dbcsr_reserve_block2d, dbcsr_scale, & dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create USE domain_submatrix_methods, ONLY: init_submatrices,& @@ -74,6 +78,7 @@ MODULE almo_scf USE kinds, ONLY: default_path_length,& dp USE mathlib, ONLY: binomial + USE message_passing, ONLY: mp_sum USE molecule_types, ONLY: get_molecule_set_info,& molecule_type USE mscfg_types, ONLY: get_matrix_from_submatrices,& @@ -116,7 +121,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_entry_scf(qs_env, calc_forces) TYPE(qs_environment_type), POINTER :: qs_env - LOGICAL :: calc_forces + LOGICAL, INTENT(IN) :: calc_forces CHARACTER(len=*), PARAMETER :: routineN = 'almo_entry_scf', routineP = moduleN//':'//routineN @@ -142,6 +147,9 @@ CONTAINS ! allow electron delocalization CALL almo_scf_delocalization(qs_env, almo_scf_env) + ! construct NLMOs + CALL construct_nlmos(qs_env, almo_scf_env) + ! electron correlation methods !CALL almo_correlation_main(qs_env,almo_scf_env) @@ -167,8 +175,8 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_init(qs_env, almo_scf_env, calc_forces) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env - LOGICAL :: calc_forces + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + LOGICAL, INTENT(IN) :: calc_forces CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_init', routineP = moduleN//':'//routineN @@ -198,6 +206,7 @@ CONTAINS almo_scf_env%opt_xalmo_diis%optimizer_type = optimizer_diis almo_scf_env%opt_xalmo_pcg%optimizer_type = optimizer_pcg almo_scf_env%opt_xalmo_trustr%optimizer_type = optimizer_trustr + almo_scf_env%opt_nlmo_pcg%optimizer_type = optimizer_pcg almo_scf_env%opt_xalmo_newton_pcg_solver%optimizer_type = optimizer_lin_eq_pcg ! get info from the qs_env @@ -508,7 +517,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_initial_guess(qs_env, almo_scf_env) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_initial_guess', & routineP = moduleN//':'//routineN @@ -760,7 +769,7 @@ CONTAINS !> \author Rustam Khaliullin ! ************************************************************************************************** SUBROUTINE almo_scf_store_extrapolation_data(almo_scf_env) - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_store_extrapolation_data', & routineP = moduleN//':'//routineN @@ -1259,7 +1268,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_main(qs_env, almo_scf_env) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_main', routineP = moduleN//':'//routineN @@ -1351,7 +1360,7 @@ CONTAINS SUBROUTINE almo_scf_delocalization(qs_env, almo_scf_env) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_delocalization', & routineP = moduleN//':'//routineN @@ -1658,6 +1667,444 @@ CONTAINS END SUBROUTINE almo_scf_delocalization +! ************************************************************************************************** +!> \brief orbital localization +!> \param qs_env ... +!> \param almo_scf_env ... +!> \par History +!> 2018.09 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE construct_nlmos(qs_env, almo_scf_env) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + + CHARACTER(len=*), PARAMETER :: routineN = 'construct_nlmos', & + routineP = moduleN//':'//routineN + + INTEGER :: ispin + + IF (almo_scf_env%construct_nlmos) THEN + + DO ispin = 1, almo_scf_env%nspins + + CALL orthogonalize_mos(ket=almo_scf_env%matrix_t(ispin), & + overlap=almo_scf_env%matrix_sigma(ispin), & + metric=almo_scf_env%matrix_s(1), & + retain_locality=.FALSE., & + only_normalize=.FALSE., & + nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), & + eps_filter=almo_scf_env%eps_filter, & + order_lanczos=almo_scf_env%order_lanczos, & + eps_lanczos=almo_scf_env%eps_lanczos, & + max_iter_lanczos=almo_scf_env%max_iter_lanczos) + ENDDO + + CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.FALSE.) + + IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%virtual_nlmos) THEN + CALL construct_virtuals(almo_scf_env) + CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.TRUE.) + ENDIF + + IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start .GT. 0.0_dp) THEN + CALL nlmo_compactification(qs_env, almo_scf_env, almo_scf_env%matrix_t) + ENDIF + + ENDIF + + END SUBROUTINE construct_nlmos + +! ************************************************************************************************** +!> \brief Calls NLMO optimization +!> \param qs_env ... +!> \param almo_scf_env ... +!> \param virtuals ... +!> \par History +!> 2019.10 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + LOGICAL, INTENT(IN) :: virtuals + + CHARACTER(len=*), PARAMETER :: routineN = 'construct_nlmos_wrapper', & + routineP = moduleN//':'//routineN + + REAL(KIND=dp) :: det_diff, prev_determinant + + almo_scf_env%overlap_determinant = 1.0 + ! KEEP: initial_vol_coeff = almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength + almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = & + -1.0_dp*almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength !NEW1 + + ! loop over the strength of the orthogonalization penalty + prev_determinant = 10.0_dp + DO WHILE (almo_scf_env%overlap_determinant .GT. almo_scf_env%opt_nlmo_pcg%opt_penalty%final_determinant) + + IF (.NOT. virtuals) THEN + CALL almo_scf_construct_nlmos(qs_env=qs_env, & + optimizer=almo_scf_env%opt_nlmo_pcg, & + matrix_s=almo_scf_env%matrix_s(1), & + matrix_mo_in=almo_scf_env%matrix_t, & + matrix_mo_out=almo_scf_env%matrix_t, & + template_matrix_sigma=almo_scf_env%matrix_sigma_inv, & + overlap_determinant=almo_scf_env%overlap_determinant, & + mat_distr_aos=almo_scf_env%mat_distr_aos, & + virtuals=virtuals, & + eps_filter=almo_scf_env%eps_filter) + ELSE + CALL almo_scf_construct_nlmos(qs_env=qs_env, & + optimizer=almo_scf_env%opt_nlmo_pcg, & + matrix_s=almo_scf_env%matrix_s(1), & + matrix_mo_in=almo_scf_env%matrix_v, & + matrix_mo_out=almo_scf_env%matrix_v, & + template_matrix_sigma=almo_scf_env%matrix_sigma_vv, & + overlap_determinant=almo_scf_env%overlap_determinant, & + mat_distr_aos=almo_scf_env%mat_distr_aos, & + virtuals=virtuals, & + eps_filter=almo_scf_env%eps_filter) + + ENDIF + + det_diff = prev_determinant - almo_scf_env%overlap_determinant + almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = & + almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength/ & + ABS(almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength_dec_factor) + + IF (det_diff < almo_scf_env%opt_nlmo_pcg%opt_penalty%determinant_tolerance) THEN + EXIT + ENDIF + prev_determinant = almo_scf_env%overlap_determinant + + ENDDO + + END SUBROUTINE construct_nlmos_wrapper + +! ************************************************************************************************** +!> \brief Construct virtual orbitals +!> \param almo_scf_env ... +!> \par History +!> 2019.10 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE construct_virtuals(almo_scf_env) + + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + + CHARACTER(len=*), PARAMETER :: routineN = 'construct_virtuals', & + routineP = moduleN//':'//routineN + + INTEGER :: ispin, n + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues + TYPE(dbcsr_type) :: tempNV1, tempVOcc1, tempVOcc2, tempVV1, & + tempVV2 + + DO ispin = 1, almo_scf_env%nspins + + CALL dbcsr_create(tempNV1, & + template=almo_scf_env%matrix_v(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempVOcc1, & + template=almo_scf_env%matrix_vo(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempVOcc2, & + template=almo_scf_env%matrix_vo(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempVV1, & + template=almo_scf_env%matrix_sigma_vv(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempVV2, & + template=almo_scf_env%matrix_sigma_vv(ispin), & + matrix_type=dbcsr_type_no_symmetry) + + ! Generate random virtual matrix + CALL dbcsr_init_random(almo_scf_env%matrix_v(ispin), & + keep_sparsity=.FALSE.) + + ! Project the orbital subspace out + CALL dbcsr_multiply("N", "N", 1.0_dp, & + almo_scf_env%matrix_s(1), & + almo_scf_env%matrix_v(ispin), & + 0.0_dp, tempNV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("T", "N", 1.0_dp, & + tempNV1, & + almo_scf_env%matrix_t(ispin), & + 0.0_dp, tempVOcc1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + tempVOcc1, & + almo_scf_env%matrix_sigma_inv(ispin), & + 0.0_dp, tempVOcc2, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("N", "T", 1.0_dp, & + almo_scf_env%matrix_t(ispin), & + tempVOcc2, & + 0.0_dp, tempNV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_add(almo_scf_env%matrix_v(ispin), tempNV1, 1.0_dp, -1.0_dp) + + ! compute VxV overlap + CALL dbcsr_multiply("N", "N", 1.0_dp, & + almo_scf_env%matrix_s(1), & + almo_scf_env%matrix_v(ispin), & + 0.0_dp, tempNV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("T", "N", 1.0_dp, & + almo_scf_env%matrix_v(ispin), & + tempNV1, & + 0.0_dp, tempVV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL orthogonalize_mos(ket=almo_scf_env%matrix_v(ispin), & + overlap=tempVV1, & + metric=almo_scf_env%matrix_s(1), & + retain_locality=.FALSE., & + only_normalize=.FALSE., & + nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), & + eps_filter=almo_scf_env%eps_filter, & + order_lanczos=almo_scf_env%order_lanczos, & + eps_lanczos=almo_scf_env%eps_lanczos, & + max_iter_lanczos=almo_scf_env%max_iter_lanczos) + + ! compute VxV block of the KS matrix + CALL dbcsr_multiply("N", "N", 1.0_dp, & + almo_scf_env%matrix_ks(ispin), & + almo_scf_env%matrix_v(ispin), & + 0.0_dp, tempNV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("T", "N", 1.0_dp, & + almo_scf_env%matrix_v(ispin), & + tempNV1, & + 0.0_dp, tempVV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_get_info(tempVV1, nfullrows_total=n) + ALLOCATE (eigenvalues(n)) + CALL cp_dbcsr_syevd(tempVV1, tempVV2, & + eigenvalues, & + para_env=almo_scf_env%para_env, & + blacs_env=almo_scf_env%blacs_env) + DEALLOCATE (eigenvalues) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + almo_scf_env%matrix_v(ispin), & + tempVV2, & + 0.0_dp, tempNV1, & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_copy(almo_scf_env%matrix_v(ispin), tempNV1) + + CALL dbcsr_release(tempNV1) + CALL dbcsr_release(tempVOcc1) + CALL dbcsr_release(tempVOcc2) + CALL dbcsr_release(tempVV1) + CALL dbcsr_release(tempVV2) + + ENDDO + + END SUBROUTINE construct_virtuals + +! ************************************************************************************************** +!> \brief Compactify (set small blocks to zero) orbitals +!> \param qs_env ... +!> \param almo_scf_env ... +!> \param matrix ... +!> \par History +!> 2019.10 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE nlmo_compactification(qs_env, almo_scf_env, matrix) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), & + INTENT(IN) :: matrix + + CHARACTER(len=*), PARAMETER :: routineN = 'nlmo_compactification', & + routineP = moduleN//':'//routineN + + INTEGER :: iblock_col, iblock_col_size, iblock_row, & + iblock_row_size, icol, irow, ispin, & + Ncols, Nrows, nspins, para_group, & + unit_nr + LOGICAL :: element_by_element + REAL(KIND=dp) :: energy, eps_local, eps_start, & + max_element, spin_factor + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ, retained + REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_iterator_type) :: iter + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_p_tmp, matrix_t_tmp + + ! define the output_unit + logger => cp_get_default_logger() + IF (logger%para_env%ionode) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + ENDIF + + nspins = SIZE(matrix) + element_by_element = .FALSE. + + IF (nspins .EQ. 1) THEN + spin_factor = 2.0_dp + ELSE + spin_factor = 1.0_dp + ENDIF + + ALLOCATE (matrix_t_tmp(nspins)) + ALLOCATE (matrix_p_tmp(nspins)) + ALLOCATE (retained(nspins)) + ALLOCATE (occ(2)) + + DO ispin = 1, nspins + + ! init temporary storage + CALL dbcsr_create(matrix_t_tmp(ispin), & + template=matrix(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_copy(matrix_t_tmp(ispin), matrix(ispin)) + + CALL dbcsr_create(matrix_p_tmp(ispin), & + template=almo_scf_env%matrix_p(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_copy(matrix_p_tmp(ispin), almo_scf_env%matrix_p(ispin)) + + ENDDO + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) + WRITE (unit_nr, '(T2,A)') & + "Energy dependence on the (block-by-block) filtering of the NLMO coefficients" + IF (unit_nr > 0) WRITE (unit_nr, '(T2,A13,A20,A20,A25)') & + "EPS filter", "Occupation Alpha", "Occupation Beta", "Energy" + ENDIF + + eps_start = almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start + eps_local = MAX(eps_start, 10E-14_dp) + + DO + + IF (eps_local > 0.11_dp) EXIT + + DO ispin = 1, nspins + + retained(ispin) = 0 + CALL dbcsr_work_create(matrix_t_tmp(ispin), work_mutable=.TRUE.) + CALL dbcsr_iterator_start(iter, matrix_t_tmp(ispin)) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, & + row_size=iblock_row_size, col_size=iblock_col_size) + DO icol = 1, iblock_col_size + + IF (element_by_element) THEN + + DO irow = 1, iblock_row_size + IF (ABS(data_p(irow, icol)) .LT. eps_local) THEN + data_p(irow, icol) = 0.0_dp + ELSE + retained(ispin) = retained(ispin) + 1 + ENDIF + ENDDO + + ELSE ! rows are blocked + + max_element = 0.0_dp + DO irow = 1, iblock_row_size + IF (ABS(data_p(irow, icol)) .GT. max_element) THEN + max_element = ABS(data_p(irow, icol)) + ENDIF + ENDDO + IF (max_element .LT. eps_local) THEN + DO irow = 1, iblock_row_size + data_p(irow, icol) = 0.0_dp + ENDDO + ELSE + retained(ispin) = retained(ispin) + iblock_row_size + ENDIF + + ENDIF ! block rows? + ENDDO ! icol + + ENDDO ! iterator + CALL dbcsr_iterator_stop(iter) + CALL dbcsr_finalize(matrix_t_tmp(ispin)) + CALL dbcsr_filter(matrix_t_tmp(ispin), eps_local) + + CALL dbcsr_get_info(matrix_t_tmp(ispin), group=para_group, & + nfullrows_total=Nrows, & + nfullcols_total=Ncols) + CALL mp_sum(retained(ispin), para_group) + + !devide by the total no. elements + occ(ispin) = retained(ispin)/Nrows/Ncols + + ! compute the global projectors (for the density matrix) + CALL almo_scf_t_to_proj( & + t=matrix_t_tmp(ispin), & + p=matrix_p_tmp(ispin), & + eps_filter=almo_scf_env%eps_filter, & + orthog_orbs=.FALSE., & + nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), & + s=almo_scf_env%matrix_s(1), & + sigma=almo_scf_env%matrix_sigma(ispin), & + sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), & + use_guess=.FALSE., & + algorithm=almo_scf_env%sigma_inv_algorithm, & + inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, & + inverse_accelerator=almo_scf_env%order_lanczos, & + eps_lanczos=almo_scf_env%eps_lanczos, & + max_iter_lanczos=almo_scf_env%max_iter_lanczos, & + para_env=almo_scf_env%para_env, & + blacs_env=almo_scf_env%blacs_env) + + ! compute dm from the projector(s) + CALL dbcsr_scale(matrix_p_tmp(ispin), spin_factor) + + ENDDO + + ! the KS matrix is updated outside the spin loop + CALL almo_dm_to_almo_ks(qs_env, & + matrix_p_tmp, & + almo_scf_env%matrix_ks, & + energy, & + almo_scf_env%eps_filter, & + almo_scf_env%mat_distr_aos) + + IF (nspins .LT. 2) occ(2) = occ(1) + IF (unit_nr > 0) WRITE (unit_nr, '(T2,E13.3,F20.10,F20.10,F25.15)') & + eps_local, occ(1), occ(2), energy + + eps_local = 2.0_dp*eps_local + + ENDDO + + DO ispin = 1, nspins + + CALL dbcsr_release(matrix_t_tmp(ispin)) + CALL dbcsr_release(matrix_p_tmp(ispin)) + + ENDDO + + DEALLOCATE (matrix_t_tmp) + DEALLOCATE (matrix_p_tmp) + DEALLOCATE (occ) + DEALLOCATE (retained) + + END SUBROUTINE nlmo_compactification + ! ***************************************************************************** !> \brief after SCF we have the final density and KS matrices compute various !> post-scf quantities @@ -1669,7 +2116,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_post(qs_env, almo_scf_env) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_post', routineP = moduleN//':'//routineN @@ -1687,6 +2134,7 @@ CONTAINS ! orthogonalize orbitals before returning them to QS ALLOCATE (matrix_t_processed(almo_scf_env%nspins)) + !ALLOCATE (matrix_v_processed(almo_scf_env%nspins)) DO ispin = 1, almo_scf_env%nspins @@ -1697,6 +2145,13 @@ CONTAINS CALL dbcsr_copy(matrix_t_processed(ispin), & almo_scf_env%matrix_t(ispin)) + !CALL dbcsr_create(matrix_v_processed(ispin), & + ! template=almo_scf_env%matrix_v(ispin), & + ! matrix_type=dbcsr_type_no_symmetry) + + !CALL dbcsr_copy(matrix_v_processed(ispin), & + ! almo_scf_env%matrix_v(ispin)) + IF (almo_scf_env%return_orthogonalized_mos) THEN CALL orthogonalize_mos(ket=matrix_t_processed(ispin), & @@ -1736,6 +2191,9 @@ CONTAINS CALL copy_dbcsr_to_fm(matrix_t_processed(ispin), mo_coeff) CALL dbcsr_release(matrix_t_processed(ispin)) ENDDO + DO ispin = 1, almo_scf_env%nspins + CALL dbcsr_release(matrix_t_processed(ispin)) + ENDDO DEALLOCATE (matrix_t_processed) ! calculate post scf properties @@ -2233,7 +2691,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_clean_up(almo_scf_env) - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_clean_up', & routineP = moduleN//':'//routineN diff --git a/src/almo_scf_env_methods.F b/src/almo_scf_env_methods.F index 56bc8b9297..e77af4e0be 100644 --- a/src/almo_scf_env_methods.F +++ b/src/almo_scf_env_methods.F @@ -15,10 +15,9 @@ MODULE almo_scf_env_methods almo_scf_env_type USE cp_control_types, ONLY: dft_control_type USE input_constants, ONLY: & - almo_constraint_distance, almo_deloc_none, almo_deloc_scf, almo_deloc_x, & - almo_deloc_x_then_scf, almo_deloc_xalmo_1diag, almo_domain_layout_atomic, & - almo_domain_layout_molecular, almo_frz_crystal, almo_mat_distr_molecular, & - almo_occ_vol_penalty_none, almo_scf_diag, almo_scf_skip, almo_scf_trustr, cg_hager_zhang, & + almo_constraint_distance, almo_deloc_none, almo_deloc_xalmo_1diag, & + almo_domain_layout_atomic, almo_domain_layout_molecular, almo_frz_crystal, & + almo_mat_distr_molecular, almo_scf_diag, almo_scf_skip, almo_scf_trustr, cg_hager_zhang, & do_bondparm_vdw, molecular_guess, tensor_orthogonal, virt_full, virt_minimal, virt_number, & xalmo_trial_r0_out USE input_section_types, ONLY: section_vals_get_subs_vals,& @@ -107,8 +106,9 @@ CONTAINS INTEGER :: handle TYPE(section_vals_type), POINTER :: almo_analysis_section, almo_opt_diis_section, & - almo_opt_pcg_section, almo_scf_section, matrix_iterate_section, penalty_section, & - xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, xalmo_opt_trustr_section + almo_opt_pcg_section, almo_scf_section, matrix_iterate_section, nlmo_opt_pcg_section, & + penalty_section, xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, & + xalmo_opt_trustr_section CALL timeset(routineN, handle) @@ -121,12 +121,13 @@ CONTAINS "XALMO_OPTIMIZER_PCG") xalmo_opt_trustr_section => section_vals_get_subs_vals(almo_scf_section, & "XALMO_OPTIMIZER_TRUSTR") + nlmo_opt_pcg_section => section_vals_get_subs_vals(almo_scf_section, & + "NLMO_OPTIMIZER_PCG") almo_analysis_section => section_vals_get_subs_vals(almo_scf_section, "ANALYSIS") xalmo_opt_newton_pcg_section => section_vals_get_subs_vals(xalmo_opt_pcg_section, & "XALMO_NEWTON_PCG_SOLVER") matrix_iterate_section => section_vals_get_subs_vals(almo_scf_section, & "MATRIX_ITERATE") - penalty_section => section_vals_get_subs_vals(almo_scf_section, "PENALTY") ! read user input ! common ALMO options @@ -154,6 +155,8 @@ CONTAINS almo_scf_env%xalmo_extrapolation_order = MAX(0, almo_scf_env%xalmo_extrapolation_order) CALL section_vals_val_get(almo_scf_section, "RETURN_ORTHOGONALIZED_MOS", & l_val=almo_scf_env%return_orthogonalized_mos) + CALL section_vals_val_get(almo_scf_section, "CONSTRUCT_NLMOS", & + l_val=almo_scf_env%construct_nlmos) CALL section_vals_val_get(matrix_iterate_section, "EPS_LANCZOS", & r_val=almo_scf_env%eps_lanczos) @@ -243,6 +246,49 @@ CONTAINS CALL section_vals_val_get(xalmo_opt_pcg_section, "PRECONDITIONER", & i_val=almo_scf_env%opt_xalmo_pcg%preconditioner) + penalty_section => section_vals_get_subs_vals(nlmo_opt_pcg_section, "PENALTY") + CALL section_vals_val_get(nlmo_opt_pcg_section, "EPS_ERROR", & + r_val=almo_scf_env%opt_nlmo_pcg%eps_error) + CALL section_vals_val_get(nlmo_opt_pcg_section, "MAX_ITER", & + i_val=almo_scf_env%opt_nlmo_pcg%max_iter) + CALL section_vals_val_get(nlmo_opt_pcg_section, "EPS_ERROR_EARLY", & + r_val=almo_scf_env%opt_nlmo_pcg%eps_error_early) + CALL section_vals_val_get(nlmo_opt_pcg_section, "MAX_ITER_EARLY", & + i_val=almo_scf_env%opt_nlmo_pcg%max_iter_early) + CALL section_vals_val_get(nlmo_opt_pcg_section, "MAX_ITER_OUTER_LOOP", & + i_val=almo_scf_env%opt_nlmo_pcg%max_iter_outer_loop) + CALL section_vals_val_get(nlmo_opt_pcg_section, "LIN_SEARCH_EPS_ERROR", & + r_val=almo_scf_env%opt_nlmo_pcg%lin_search_eps_error) + CALL section_vals_val_get(nlmo_opt_pcg_section, "LIN_SEARCH_STEP_SIZE_GUESS", & + r_val=almo_scf_env%opt_nlmo_pcg%lin_search_step_size_guess) + CALL section_vals_val_get(nlmo_opt_pcg_section, "PRECOND_FILTER_THRESHOLD", & + r_val=almo_scf_env%opt_nlmo_pcg%neglect_threshold) + CALL section_vals_val_get(nlmo_opt_pcg_section, "CONJUGATOR", & + i_val=almo_scf_env%opt_nlmo_pcg%conjugator) + CALL section_vals_val_get(nlmo_opt_pcg_section, "PRECONDITIONER", & + i_val=almo_scf_env%opt_nlmo_pcg%preconditioner) + CALL section_vals_val_get(penalty_section, & + "OPERATOR", & + i_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%operator_type) + CALL section_vals_val_get(penalty_section, & + "PENALTY_STRENGTH", & + r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength) + CALL section_vals_val_get(penalty_section, & + "PENALTY_STRENGTH_DECREASE_FACTOR", & + r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength_dec_factor) + CALL section_vals_val_get(penalty_section, & + "DETERMINANT_TOLERANCE", & + r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%determinant_tolerance) + CALL section_vals_val_get(penalty_section, & + "FINAL_DETERMINANT", & + r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%final_determinant) + CALL section_vals_val_get(penalty_section, & + "COMPACTIFICATION_FILTER_START", & + r_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start) + CALL section_vals_val_get(penalty_section, & + "VIRTUAL_NLMOS", & + l_val=almo_scf_env%opt_nlmo_pcg%opt_penalty%virtual_nlmos) + CALL section_vals_val_get(xalmo_opt_newton_pcg_section, "EPS_ERROR", & r_val=almo_scf_env%opt_xalmo_newton_pcg_solver%eps_error) CALL section_vals_val_get(xalmo_opt_newton_pcg_section, "MAX_ITER", & @@ -252,13 +298,6 @@ CONTAINS CALL section_vals_val_get(xalmo_opt_newton_pcg_section, "PRECONDITIONER", & i_val=almo_scf_env%opt_xalmo_newton_pcg_solver%preconditioner) - CALL section_vals_val_get(penalty_section, & - "OCCUPIED_VOLUME_PENALTY_COEFF", & - r_val=almo_scf_env%penalty%occ_vol_coeff) - CALL section_vals_val_get(penalty_section, & - "OCCUPIED_VOLUME_PENALTY_METHOD", & - i_val=almo_scf_env%penalty%occ_vol_method) - CALL section_vals_val_get(almo_analysis_section, "_SECTION_PARAMETERS_", & l_val=almo_scf_env%almo_analysis%do_analysis) CALL section_vals_val_get(almo_analysis_section, "FROZEN_MO_ENERGY_TERM", & @@ -504,16 +543,6 @@ CONTAINS ENDIF ! end analysis settings - ! check penalty settings - IF (almo_scf_env%penalty%occ_vol_method .NE. almo_occ_vol_penalty_none) THEN - IF (almo_scf_env%deloc_method .NE. almo_deloc_x .AND. & - almo_scf_env%deloc_method .NE. almo_deloc_scf .AND. & - almo_scf_env%deloc_method .NE. almo_deloc_x_then_scf) THEN - CALL cp_abort(__LOCATION__, & - "Occupied volume penalty seems to work only with completely delocalized orbitals") - ENDIF - ENDIF ! end penalty settings - CALL timestop(handle) END SUBROUTINE almo_scf_init_read_write_input diff --git a/src/almo_scf_lbfgs_types.F b/src/almo_scf_lbfgs_types.F new file mode 100644 index 0000000000..139fc6b266 --- /dev/null +++ b/src/almo_scf_lbfgs_types.F @@ -0,0 +1,330 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright (C) 2000 - 2020 CP2K developers group ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Limited memory BFGS +!> \par History +!> 2019.10 created [Rustam Z Khaliullin] +!> \author Rustam Z Khaliullin +! ************************************************************************************************** +MODULE almo_scf_lbfgs_types + !USE cp_external_control, ONLY: external_control + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_unit_nr + !USE cp_log_handling, ONLY: cp_to_string + USE dbcsr_api, ONLY: dbcsr_add,& + dbcsr_copy,& + dbcsr_create,& + dbcsr_dot,& + dbcsr_release,& + dbcsr_scale,& + dbcsr_type + USE kinds, ONLY: dp +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_lbfgs_types' + + PUBLIC :: lbfgs_seed, & + lbfgs_create, & + lbfgs_release, & + lbfgs_get_direction, & + lbfgs_history_type + + TYPE lbfgs_history_type + INTEGER :: nstore + ! istore counts the total number of action=2 pushes + ! istore is designed to become more than nstore eventually + ! there are two counters: the main variable and gradient + INTEGER, DIMENSION(2) :: istore + TYPE(dbcsr_type), DIMENSION(:, :, :), ALLOCATABLE :: matrix + REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: rho + END TYPE lbfgs_history_type + +CONTAINS + +! ************************************************************************************************** +!> \brief interface subroutine to store the first variable/gradient pair +!> \param history ... +!> \param variable ... +!> \param gradient ... +! ************************************************************************************************** + SUBROUTINE lbfgs_seed(history, variable, gradient) + + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: variable, gradient + + CALL lbfgs_history_push(history, variable, vartype=1, action=1) + CALL lbfgs_history_push(history, gradient, vartype=2, action=1) + + END SUBROUTINE lbfgs_seed + +! ************************************************************************************************** +!> \brief interface subroutine to store a variable/gradient pair +!> and predict direction +!> \param history ... +!> \param variable ... +!> \param gradient ... +!> \param direction ... +! ************************************************************************************************** + + SUBROUTINE lbfgs_get_direction(history, variable, gradient, direction) + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: variable, gradient + TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: direction + + ! action 2 will calculate delta = (new - old) + ! in the last used storage cell + CALL lbfgs_history_push(history, variable, vartype=1, action=2) + CALL lbfgs_history_push(history, gradient, vartype=2, action=2) + ! compute rho for the last stored value + CALL lbfgs_history_last_rho(history) + + CALL lbfgs_history_direction(history, gradient, direction) + + ! action 1 will seed the next storage cell + CALL lbfgs_history_push(history, variable, vartype=1, action=1) + CALL lbfgs_history_push(history, gradient, vartype=2, action=1) + + END SUBROUTINE lbfgs_get_direction + +! ************************************************************************************************** +!> \brief create history storage for limited memory bfgs +!> \param history ... +!> \param nspins ... +!> \param nstore ... +! ************************************************************************************************** + SUBROUTINE lbfgs_create(history, nspins, nstore) + + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + INTEGER, INTENT(IN) :: nspins, nstore + + INTEGER :: nallocate + + nallocate = MAX(1, nstore) + history%nstore = nallocate + history%istore(:) = 0 ! total number of action-2 pushes + ALLOCATE (history%matrix(nspins, nallocate, 2)) + ALLOCATE (history%rho(nspins, nallocate)) + + END SUBROUTINE lbfgs_create + +! ************************************************************************************************** +!> \brief release the bfgs history +!> \param history ... +! ************************************************************************************************** + SUBROUTINE lbfgs_release(history) + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + + INTEGER :: ispin, istore, ivartype + + ! delete history + DO ispin = 1, SIZE(history%matrix, 1) + DO ivartype = 1, 2 + DO istore = 1, MIN(history%istore(ivartype) + 1, history%nstore) + !WRITE(*,*) "ZREL: ispin,istore,vartype", ispin, istore, ivartype + CALL dbcsr_release(history%matrix(ispin, istore, ivartype)) + ENDDO + ENDDO + ENDDO + DEALLOCATE (history%matrix) + DEALLOCATE (history%rho) + + END SUBROUTINE lbfgs_release + +! ************************************************************************************************** +!> \brief once all data in the last cell is stored, compute rho +!> \param history ... +! ************************************************************************************************** + SUBROUTINE lbfgs_history_last_rho(history) + + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + + INTEGER :: ispin, istore + + !logger => cp_get_default_logger() + !IF (logger%para_env%mepos == logger%para_env%source) THEN + ! unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + !ELSE + ! unit_nr = -1 + !ENDIF + + DO ispin = 1, SIZE(history%matrix, 1) + + istore = MOD(history%istore(1) - 1, history%nstore) + 1 + CALL dbcsr_dot(history%matrix(ispin, istore, 1), & + history%matrix(ispin, istore, 2), & + history%rho(ispin, istore)) + + history%rho(ispin, istore) = 1.0_dp/history%rho(ispin, istore) + + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, *) "Rho in cell ", istore, " is computed ", history%rho(ispin, istore) + !ENDIF + + ENDDO ! ispin + + END SUBROUTINE lbfgs_history_last_rho + +! ************************************************************************************************** +!> \brief store data in history +!> vartype - which data piece to store: 1 - variable, 2 - gradient +!> operation - what to do: 1 - erase existing and store new +!> 2 - store = new - existing +!> \param history ... +!> \param matrix ... +!> \param vartype ... +!> \param action ... +! ************************************************************************************************** + SUBROUTINE lbfgs_history_push(history, matrix, vartype, action) + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: matrix + INTEGER, INTENT(IN) :: vartype, action + + INTEGER :: ispin, istore + + !logger => cp_get_default_logger() + !IF (logger%para_env%mepos == logger%para_env%source) THEN + ! unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + !ELSE + ! unit_nr = -1 + !ENDIF + + ! increase the counter: it moves the pointer to the next cell + ! for action==1 this is a "pretend" increase; the pointer will be moved back in the end + history%istore(vartype) = history%istore(vartype) + 1 + + DO ispin = 1, SIZE(history%matrix, 1) + + istore = MOD(history%istore(vartype) - 1, history%nstore) + 1 + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, *) "Action ", action, " modifying cell ", istore + !END IF + + IF (history%istore(vartype) <= history%nstore .AND. & + action .EQ. 1) THEN + !WRITE(*,*) "ZCRE: ispin,istore,vartype", ispin, istore, vartype + CALL dbcsr_create(history%matrix(ispin, istore, vartype), & + template=matrix(ispin)) + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, *) "Creating new matrix..." + !END IF + ENDIF + + IF (action .EQ. 1) THEN + CALL dbcsr_copy(history%matrix(ispin, istore, vartype), matrix(ispin)) + ELSE + CALL dbcsr_add(history%matrix(ispin, istore, vartype), matrix(ispin), -1.0_dp, 1.0_dp) + ENDIF + + ENDDO ! ispin + + ! allow the pointer to move forward only if deltas are stored (action==2) + ! otherwise return the pointer to the previous value + IF (action .EQ. 1) THEN + history%istore(vartype) = history%istore(vartype) - 1 + ENDIF + + END SUBROUTINE lbfgs_history_push + +! ************************************************************************************************** +!> \brief use history data to construct dir = -Hinv.grad +!> \param history ... +!> \param gradient ... +!> \param direction ... +! ************************************************************************************************** + SUBROUTINE lbfgs_history_direction(history, gradient, direction) + + TYPE(lbfgs_history_type), INTENT(INOUT) :: history + TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: gradient + TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: direction + + INTEGER :: ispin, istore, iterm, nterms + REAL(KIND=dp) :: beta, gammak + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha + TYPE(dbcsr_type) :: q + + !logger => cp_get_default_logger() + !IF (logger%para_env%mepos == logger%para_env%source) THEN + ! unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + !ELSE + ! unit_nr = -1 + !ENDIF + + IF (history%istore(1) .NE. history%istore(2)) THEN + CPABORT("BFGS APIs are not used correctly") + ENDIF + + nterms = MIN(history%istore(1), history%nstore) + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, *) "L-BFGS terms used: ", nterms + !END IF + + ALLOCATE (alpha(nterms)) + + DO ispin = 1, SIZE(history%matrix, 1) + + CALL dbcsr_create(q, template=gradient(ispin)) + + CALL dbcsr_copy(q, gradient(ispin)) + + ! loop over all stored items + DO iterm = 1, nterms + + ! location: from recent to oldest stored + istore = MOD(history%istore(1) - iterm, history%nstore) + 1 + + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, *) "Record locator: ", istore + !END IF + + CALL dbcsr_dot(history%matrix(ispin, istore, 1), q, alpha(iterm)) + alpha(iterm) = history%rho(ispin, istore)*alpha(iterm) + CALL dbcsr_add(q, history%matrix(ispin, istore, 2), 1.0_dp, -alpha(iterm)) + + ! use the most recent term to + ! compute gamma_k, Nocedal (7.20) and then get H0 + IF (iterm .EQ. 1) THEN + CALL dbcsr_dot(history%matrix(ispin, istore, 2), history%matrix(ispin, istore, 2), gammak) + gammak = 1.0_dp/(gammak*history%rho(ispin, istore)) + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, *) "Gamma_k: ", gammak + !END IF + ENDIF + + ENDDO ! iterm, first loop from recent to oldest + + ! now q stores Nocedal's r = (gamma_k*I).q + CALL dbcsr_scale(q, gammak) + + ! loop over all stored items + DO iterm = nterms, 1, -1 + + ! location: from oldest to recent stored + istore = MOD(history%istore(1) - iterm, history%nstore) + 1 + + CALL dbcsr_dot(history%matrix(ispin, istore, 2), q, beta) + beta = history%rho(ispin, istore)*beta + CALL dbcsr_add(q, history%matrix(ispin, istore, 1), 1.0_dp, alpha(iterm) - beta) + + ENDDO ! iterm, forst loop from recent to oldest + + !RZK-warning: unclear whether q should be multiplied by minus one + CALL dbcsr_scale(q, -1.0) + CALL dbcsr_copy(direction(ispin), q) + + CALL dbcsr_release(q) + + ENDDO !ispin + + DEALLOCATE (alpha) + + END SUBROUTINE lbfgs_history_direction + +END MODULE almo_scf_lbfgs_types + diff --git a/src/almo_scf_methods.F b/src/almo_scf_methods.F index c478003414..c6c0d4399b 100644 --- a/src/almo_scf_methods.F +++ b/src/almo_scf_methods.F @@ -75,10 +75,53 @@ MODULE almo_scf_methods distribute_domains, & almo_scf_ks_to_ks_xx, & construct_domain_r_down, & - xalmo_initial_guess + xalmo_initial_guess, & + fill_matrix_with_ones CONTAINS +! ************************************************************************************************** +!> \brief Fill all matrix blocks with 1.0_dp +!> \param matrix ... +!> \par History +!> 2019.09 created [Rustam Z Khaliullin] +!> \author Rustam Z Khaliullin +! ************************************************************************************************** + SUBROUTINE fill_matrix_with_ones(matrix) + + TYPE(dbcsr_type), INTENT(INOUT) :: matrix + + INTEGER :: col, hold, iblock_col, iblock_row, & + mynode, nblkcols_tot, nblkrows_tot, row + LOGICAL :: tr + REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_new_block + TYPE(dbcsr_distribution_type) :: dist + + CALL dbcsr_get_info(matrix, distribution=dist) + CALL dbcsr_distribution_get(dist, mynode=mynode) + CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) + nblkrows_tot = dbcsr_nblkrows_total(matrix) + nblkcols_tot = dbcsr_nblkcols_total(matrix) + DO row = 1, nblkrows_tot + DO col = 1, nblkcols_tot + tr = .FALSE. + iblock_row = row + iblock_col = col + CALL dbcsr_get_stored_coordinates(matrix, & + iblock_row, iblock_col, hold) + IF (hold .EQ. mynode) THEN + NULLIFY (p_new_block) + CALL dbcsr_reserve_block2d(matrix, & + iblock_row, iblock_col, p_new_block) + CPASSERT(ASSOCIATED(p_new_block)) + p_new_block(:, :) = 1.0_dp + ENDIF + ENDDO + ENDDO + CALL dbcsr_finalize(matrix) + + END SUBROUTINE fill_matrix_with_ones + ! ************************************************************************************************** !> \brief builds projected KS matrices for the overlapping domains !> also computes the DIIS error vector as a by-product @@ -2198,7 +2241,7 @@ CONTAINS TYPE(domain_submatrix_type), DIMENSION(:), & INTENT(IN), OPTIONAL :: subm_s_inv, subm_s_inv_half, & subm_s_half, subm_r_down - TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: matrix_trimmer + TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_trimmer TYPE(dbcsr_type), INTENT(IN) :: dpattern TYPE(domain_map_type), INTENT(IN) :: map INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain diff --git a/src/almo_scf_optimizer.F b/src/almo_scf_optimizer.F index ea64e67a8b..eae7630fd5 100644 --- a/src/almo_scf_optimizer.F +++ b/src/almo_scf_optimizer.F @@ -16,18 +16,25 @@ MODULE almo_scf_optimizer almo_scf_diis_push,& almo_scf_diis_release,& almo_scf_diis_type + USE almo_scf_lbfgs_types, ONLY: lbfgs_create,& + lbfgs_get_direction,& + lbfgs_history_type,& + lbfgs_release,& + lbfgs_seed USE almo_scf_methods, ONLY: & almo_scf_ks_blk_to_tv_blk, almo_scf_ks_to_ks_blk, almo_scf_ks_to_ks_xx, & almo_scf_ks_xx_to_tv_xx, almo_scf_p_blk_to_t_blk, almo_scf_t_rescaling, & almo_scf_t_to_proj, apply_domain_operators, apply_projector, & construct_domain_preconditioner, construct_domain_r_down, construct_domain_s_inv, & - construct_domain_s_sqrt, get_overlap, orthogonalize_mos, pseudo_invert_diagonal_blk, & - xalmo_initial_guess + construct_domain_s_sqrt, fill_matrix_with_ones, get_overlap, orthogonalize_mos, & + pseudo_invert_diagonal_blk, xalmo_initial_guess USE almo_scf_qs, ONLY: almo_dm_to_almo_ks,& almo_dm_to_qs_env,& - almo_scf_update_ks_energy + almo_scf_update_ks_energy,& + matrix_qs_to_almo USE almo_scf_types, ONLY: almo_scf_env_type,& optimizer_options_type + USE cell_types, ONLY: cell_type USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& cp_dbcsr_cholesky_invert,& @@ -37,7 +44,8 @@ MODULE almo_scf_optimizer open_file USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_unit_nr,& - cp_logger_type + cp_logger_type,& + cp_to_string USE cp_output_handling, ONLY: cp_print_key_finished_output,& cp_print_key_unit_nr USE cp_para_types, ONLY: cp_para_env_type @@ -70,11 +78,11 @@ MODULE almo_scf_optimizer domain_submatrix_type,& select_row USE input_constants, ONLY: & - almo_occ_vol_penalty_none, almo_scf_diag, almo_scf_dm_sign, cg_dai_yuan, cg_fletcher, & - cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, cg_liu_storey, cg_polak_ribiere, & - cg_zero, trustr_cauchy, trustr_dogleg, virt_full, xalmo_case_block_diag, & - xalmo_case_fully_deloc, xalmo_case_normal, xalmo_prec_domain, xalmo_prec_full, & - xalmo_prec_zero + almo_scf_diag, almo_scf_dm_sign, cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, & + cg_hager_zhang, cg_hestenes_stiefel, cg_liu_storey, cg_polak_ribiere, cg_zero, & + op_loc_berry, op_loc_pipek, trustr_cauchy, trustr_dogleg, virt_full, & + xalmo_case_block_diag, xalmo_case_fully_deloc, xalmo_case_normal, xalmo_prec_domain, & + xalmo_prec_full, xalmo_prec_zero USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type USE iterate_matrix, ONLY: determinant,& @@ -83,9 +91,15 @@ MODULE almo_scf_optimizer USE kinds, ONLY: dp USE machine, ONLY: m_flush,& m_walltime + USE message_passing, ONLY: mp_sum + USE particle_methods, ONLY: get_particle_set + USE particle_types, ONLY: particle_type USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type + USE qs_kind_types, ONLY: qs_kind_type + USE qs_loc_utils, ONLY: compute_berry_operator + USE qs_localization_methods, ONLY: initialize_weights #include "./base/base_uses.f90" IMPLICIT NONE @@ -97,7 +111,8 @@ MODULE almo_scf_optimizer PUBLIC :: almo_scf_block_diagonal, & almo_scf_xalmo_eigensolver, & almo_scf_xalmo_trustr, & - almo_scf_xalmo_pcg + almo_scf_xalmo_pcg, & + almo_scf_construct_nlmos LOGICAL, PARAMETER :: debug_mode = .FALSE. LOGICAL, PARAMETER :: safe_mode = .FALSE. @@ -119,8 +134,8 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_block_diagonal(qs_env, almo_scf_env, optimizer) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env - TYPE(optimizer_options_type) :: optimizer + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + TYPE(optimizer_options_type), INTENT(IN) :: optimizer CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_block_diagonal', & routineP = moduleN//':'//routineN @@ -138,8 +153,6 @@ CONTAINS TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_mixing_old_blk TYPE(qs_energy_type), POINTER :: qs_energy -!TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks - CALL timeset(routineN, handle) ! get a useful output_unit @@ -438,8 +451,8 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, optimizer) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env - TYPE(optimizer_options_type) :: optimizer + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + TYPE(optimizer_options_type), INTENT(IN) :: optimizer CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_xalmo_eigensolver', & routineP = moduleN//':'//routineN @@ -855,9 +868,10 @@ CONTAINS special_case) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env - TYPE(optimizer_options_type) :: optimizer - TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: quench_t, matrix_t_in, matrix_t_out + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + TYPE(optimizer_options_type), INTENT(IN) :: optimizer + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), & + INTENT(INOUT) :: quench_t, matrix_t_in, matrix_t_out LOGICAL, INTENT(IN) :: assume_t0_q0x, perturbation_only INTEGER, INTENT(IN), OPTIONAL :: special_case @@ -865,24 +879,28 @@ CONTAINS routineP = moduleN//':'//routineN CHARACTER(LEN=20) :: iter_type - INTEGER :: cg_iteration, fixed_line_search_niter, handle, ispin, iteration, & - line_search_iteration, max_iter, my_special_case, ndomains, nspins, outer_iteration, & - outer_max_iter, prec_type, unit_nr + INTEGER :: cg_iteration, dim_op, fixed_line_search_niter, handle, idim0, ielem, ispin, & + iteration, line_search_iteration, max_iter, my_special_case, ndomains, nmo, nspins, & + outer_iteration, outer_max_iter, para_group, prec_type, reim, unit_nr INTEGER, ALLOCATABLE, DIMENSION(:) :: nocc LOGICAL :: blissful_neglect, converged, just_started, line_search, normalize_orbitals, & - optimize_theta, outer_prepare_to_exit, penalty_occ_vol, prepare_to_exit, & - reset_conjugator, skip_grad, use_guess - REAL(kind=dp) :: appr_sec_der, beta, denom, denom2, e0, e1, energy_diff, energy_new, & - energy_old, eps_skip_gradients, g0, g1, grad_norm, grad_norm_frob, line_search_error, & - next_step_size_guess, penalty_amplitude, penalty_func_new, spin_factor, step_size, t1, & - t2, tempreal + optimize_theta, outer_prepare_to_exit, penalty_occ_local, penalty_occ_vol, & + prepare_to_exit, reset_conjugator, skip_grad, use_guess + REAL(dp), ALLOCATABLE, DIMENSION(:) :: reim_diag, weights, z2 + REAL(kind=dp) :: appr_sec_der, beta, denom, denom2, e0, e1, energy_coeff, energy_diff, & + energy_new, energy_old, eps_skip_gradients, fval, g0, g1, grad_norm, grad_norm_frob, & + line_search_error, localiz_coeff, localization_obj_function, next_step_size_guess, & + penalty_amplitude, penalty_func_new, spin_factor, step_size, t1, t2, tempreal REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: grad_norm_spin, & penalty_occ_vol_g_prefactor, & penalty_occ_vol_h_prefactor + TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: qs_matrix_s + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set_almo, op_sm_set_qs TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: FTsiginv, grad, m_sig_sqrti_ii, m_t_in_local, & m_theta, prec_vv, prev_grad, prev_minus_prec_grad, prev_step, siginvTFTsiginv, ST, step, & - STsiginv_0 + STsiginv_0, tempNOcc, tempNOcc_1, tempOccOcc TYPE(domain_submatrix_type), ALLOCATABLE, & DIMENSION(:, :) :: bad_modes_projector_down, domain_r_down @@ -938,15 +956,14 @@ CONTAINS eps_skip_gradients = almo_scf_env%real01 ! penalty amplitude adjusts the strenght of volume conservation - ! the following guidelines are useful - ! A = T for n = 2 - ! A = 2T for n = 4 - ! A = (32/6)T for n = 6 - penalty_occ_vol = (almo_scf_env%penalty%occ_vol_method .NE. almo_occ_vol_penalty_none .AND. & - my_special_case .EQ. xalmo_case_fully_deloc) - normalize_orbitals = penalty_occ_vol - ! not used with lndet: penalty_order=2 - penalty_amplitude = almo_scf_env%penalty%occ_vol_coeff + energy_coeff = 1.0_dp !optimizer%opt_penalty%energy_coeff + localiz_coeff = 0.0_dp !optimizer%opt_penalty%occ_loc_coeff + penalty_amplitude = 0.0_dp !optimizer%opt_penalty%occ_vol_coeff + penalty_occ_vol = .FALSE. !( optimizer%opt_penalty%occ_vol_method & + !.NE. penalty_type_none .AND. my_special_case .EQ. xalmo_case_fully_deloc ) + penalty_occ_local = .FALSE. !( optimizer%opt_penalty%occ_loc_method & + !.NE. penalty_type_none .AND. my_special_case .EQ. xalmo_case_fully_deloc ) + normalize_orbitals = penalty_occ_vol .OR. penalty_occ_local ALLOCATE (penalty_occ_vol_g_prefactor(nspins)) ALLOCATE (penalty_occ_vol_h_prefactor(nspins)) penalty_occ_vol_g_prefactor(:) = 0.0_dp @@ -989,6 +1006,48 @@ CONTAINS matrix_type=dbcsr_type_no_symmetry) ENDDO + ! Compute localization matrices + IF (penalty_occ_local) THEN + + CALL get_qs_env(qs_env=qs_env, & + matrix_s=qs_matrix_s, & + cell=cell) + + IF (cell%orthorhombic) THEN + dim_op = 3 + ELSE + dim_op = 6 + END IF + ALLOCATE (weights(6)) + weights = 0.0_dp + + CALL initialize_weights(cell, weights) + + ALLOCATE (op_sm_set_qs(2, dim_op)) + ALLOCATE (op_sm_set_almo(2, dim_op)) + + DO idim0 = 1, dim_op + DO reim = 1, SIZE(op_sm_set_qs, 1) + NULLIFY (op_sm_set_qs(reim, idim0)%matrix) + ALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, & + name="almo_scf_env%op_sm_"//TRIM(ADJUSTL(cp_to_string(reim)))//"-"//TRIM(ADJUSTL(cp_to_string(idim0)))) + CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp) + NULLIFY (op_sm_set_almo(reim, idim0)%matrix) + ALLOCATE (op_sm_set_almo(reim, idim0)%matrix) + CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%matrix_s(1), & + name="almo_scf_env%op_sm_"//TRIM(ADJUSTL(cp_to_string(reim)))//"-"//TRIM(ADJUSTL(cp_to_string(idim0)))) + CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp) + ENDDO + END DO + + CALL compute_berry_operator(qs_env, cell, op_sm_set_qs, dim_op) + + !CALL matrix_qs_to_almo(op_sm_set_qs, op_sm_set_almo, & + ! almo_scf_env%mat_distr_aos, .FALSE.) + + ENDIF + ! create initial guess from the initial orbitals CALL xalmo_initial_guess(m_guess=m_theta, & m_t_in=m_t_in_local, & @@ -1024,6 +1083,9 @@ CONTAINS ALLOCATE (step(nspins)) ALLOCATE (prev_minus_prec_grad(nspins)) ALLOCATE (m_sig_sqrti_ii(nspins)) + ALLOCATE (tempNOcc(nspins)) + ALLOCATE (tempNOcc_1(nspins)) + ALLOCATE (tempOccOcc(nspins)) DO ispin = 1, nspins ! init temporary storage @@ -1060,6 +1122,15 @@ CONTAINS CALL dbcsr_create(m_sig_sqrti_ii(ispin), & template=almo_scf_env%matrix_sigma_inv(ispin), & matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempNOcc(ispin), & + template=matrix_t_out(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempNOcc_1(ispin), & + template=matrix_t_out(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc(ispin), & + template=almo_scf_env%matrix_sigma_inv(ispin), & + matrix_type=dbcsr_type_no_symmetry) CALL dbcsr_set(step(ispin), 0.0_dp) CALL dbcsr_set(prev_step(ispin), 0.0_dp) @@ -1119,6 +1190,56 @@ CONTAINS ENDIF ! assume_t0_q0x + ! localization functional + IF (penalty_occ_local) THEN + + ! compute S.R0.B.R0.S + CALL dbcsr_multiply("N", "N", 1.0_dp, & + almo_scf_env%matrix_s(1), & + matrix_t_in(ispin), & + 0.0_dp, tempNOcc(ispin), & + filter_eps=almo_scf_env%eps_filter) + CALL dbcsr_multiply("N", "N", 1.0_dp, & + tempNOcc(ispin), & + almo_scf_env%matrix_sigma_inv(ispin), & + 0.0_dp, tempNOCC_1(ispin), & + filter_eps=almo_scf_env%eps_filter) + + DO idim0 = 1, SIZE(op_sm_set_qs, 2) ! this loop is over miller ind + DO reim = 1, SIZE(op_sm_set_qs, 1) ! this loop is over Re/Im + + CALL matrix_qs_to_almo(op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix, & + almo_scf_env%mat_distr_aos, .FALSE.) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + op_sm_set_almo(reim, idim0)%matrix, & + matrix_t_in(ispin), & + 0.0_dp, tempNOcc(ispin), & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("T", "N", 1.0_dp, & + matrix_t_in(ispin), & + tempNOcc(ispin), & + 0.0_dp, tempOccOcc(ispin), & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + tempNOCC_1(ispin), & + tempOccOcc(ispin), & + 0.0_dp, tempNOcc(ispin), & + filter_eps=almo_scf_env%eps_filter) + + CALL dbcsr_multiply("N", "T", 1.0_dp, & + tempNOcc(ispin), & + tempNOcc_1(ispin), & + 0.0_dp, op_sm_set_almo(reim, idim0)%matrix, & + filter_eps=almo_scf_env%eps_filter) + + ENDDO + ENDDO ! end loop over idim0 + + ENDIF !penalty_occ_local + ENDDO ! ispin ! start the outer SCF loop @@ -1141,6 +1262,8 @@ CONTAINS line_search_iteration = 0 energy_new = 0.0_dp energy_old = 0.0_dp + energy_diff = 0.0_dp + localization_obj_function = 0.0_dp line_search_error = 0.0_dp t1 = m_walltime() @@ -1180,7 +1303,69 @@ CONTAINS -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin) penalty_occ_vol_h_prefactor(ispin) = 0.0_dp ENDIF - ENDDO ! ispin + ENDDO + + localization_obj_function = 0.0_dp + ! RZK-warning: This block must be combined with the loss function + IF (penalty_occ_local) THEN + DO ispin = 1, nspins + + ! LzL insert localization penalty + localization_obj_function = 0.0_dp + CALL dbcsr_get_info(almo_scf_env%matrix_sigma_inv(ispin), nfullrows_total=nmo) + ALLOCATE (z2(nmo)) + ALLOCATE (reim_diag(nmo)) + + CALL dbcsr_get_info(tempOccOcc(ispin), group=para_group) + + DO idim0 = 1, SIZE(op_sm_set_qs, 2) ! this loop is over miller ind + + z2(:) = 0.0_dp + + DO reim = 1, SIZE(op_sm_set_qs, 1) ! this loop is over Re/Im + + !CALL matrix_qs_to_almo(op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix, & + ! almo_scf_env%mat_distr_aos, .FALSE.) + CALL dbcsr_multiply("N", "N", 1.0_dp, & + op_sm_set_almo(reim, idim0)%matrix, & + matrix_t_out(ispin), & + 0.0_dp, tempNOcc(ispin), & + filter_eps=almo_scf_env%eps_filter) + !warning - save time by computing only the diagonal elements + CALL dbcsr_multiply("T", "N", 1.0_dp, & + matrix_t_out(ispin), & + tempNOcc(ispin), & + 0.0_dp, tempOccOcc(ispin), & + filter_eps=almo_scf_env%eps_filter) + + reim_diag = 0.0_dp + CALL dbcsr_get_diag(tempOccOcc(ispin), reim_diag) + CALL mp_sum(reim_diag, para_group) + z2(:) = z2(:) + reim_diag(:)*reim_diag(:) + + ENDDO + + DO ielem = 1, nmo + SELECT CASE (2) ! allows for selection of different spread functionals + CASE (1) ! functional = -W_I * log( |z_I|^2 ) + fval = -weights(idim0)*LOG(ABS(z2(ielem))) + CASE (2) ! functional = W_I * ( 1 - |z_I|^2 ) + fval = weights(idim0) - weights(idim0)*ABS(z2(ielem)) + CASE (3) ! functional = W_I * ( 1 - |z_I| ) + fval = weights(idim0) - weights(idim0)*SQRT(ABS(z2(ielem))) + END SELECT + localization_obj_function = localization_obj_function + fval + ENDDO + + ENDDO ! end loop over idim0 + + DEALLOCATE (z2) + DEALLOCATE (reim_diag) + + energy_new = energy_new + localiz_coeff*localization_obj_function + + ENDDO ! ispin + ENDIF ! penalty_occ_local DO ispin = 1, nspins @@ -1209,6 +1394,7 @@ CONTAINS DO ispin = 1, nspins + ! RZK-warning: z2 supplied as input seems to be deallocated CALL compute_gradient( & m_grad_out=grad(ispin), & m_ks=almo_scf_env%matrix_ks(ispin), & @@ -1235,7 +1421,13 @@ CONTAINS envelope_amplitude=almo_scf_env%envelope_amplitude, & eps_filter=almo_scf_env%eps_filter, & spin_factor=spin_factor, & - special_case=my_special_case) + special_case=my_special_case, & + penalty_occ_local=penalty_occ_local, & + op_sm_set=op_sm_set_almo, & + weights=weights, & + energy_coeff=energy_coeff, & + localiz_coeff=localiz_coeff, & + z2=z2) ENDDO ! ispin @@ -1298,9 +1490,6 @@ CONTAINS DO ispin = 1, nspins CALL dbcsr_norm(grad(ispin), dbcsr_norm_maxabsnorm, & norm_scalar=grad_norm_spin(ispin)) - !grad_norm_frob = dbcsr_frobenius_norm(grad(ispin)) / & - ! dbcsr_frobenius_norm(quench_t(ispin)) - !IF (unit_nr > 0 ) WRITE(*,*) "Gradient RMS norm: ", grad_norm_frob ENDDO ! ispin grad_norm = MAXVAL(grad_norm_spin) @@ -1674,19 +1863,25 @@ CONTAINS iter_type, iteration, & energy_new, energy_diff, grad_norm, & t2 - t1 + IF (penalty_occ_local .OR. penalty_occ_vol) THEN + WRITE (unit_nr, '(T2,A25,F23.10)') & + "Energy component:", (energy_new - penalty_func_new - localization_obj_function) + ENDIF + IF (penalty_occ_local) THEN + WRITE (unit_nr, '(T2,A25,F23.10)') & + "Localization component:", localization_obj_function + ENDIF IF (penalty_occ_vol) THEN - WRITE (unit_nr, '(T2,A19,F23.10)') & - "Energy component:", energy_new - penalty_func_new - WRITE (unit_nr, '(T2,A19,F23.10)') & + WRITE (unit_nr, '(T2,A25,F23.10)') & "Penalty component:", penalty_func_new ENDIF ENDIF IF (my_special_case .EQ. xalmo_case_block_diag) THEN IF (penalty_occ_vol) THEN - almo_scf_env%almo_scf_energy = energy_new - penalty_func_new + almo_scf_env%almo_scf_energy = energy_new - penalty_func_new - localization_obj_function ELSE - almo_scf_env%almo_scf_energy = energy_new + almo_scf_env%almo_scf_energy = energy_new - localization_obj_function ENDIF ENDIF @@ -1742,8 +1937,14 @@ CONTAINS CALL dbcsr_release(m_sig_sqrti_ii(ispin)) CALL release_submatrices(domain_r_down(:, ispin)) CALL release_submatrices(bad_modes_projector_down(:, ispin)) + CALL dbcsr_release(tempNOcc(ispin)) + CALL dbcsr_release(tempNOcc_1(ispin)) + CALL dbcsr_release(tempOccOcc(ispin)) ENDDO ! ispin + DEALLOCATE (tempNOcc) + DEALLOCATE (tempNOcc_1) + DEALLOCATE (tempOccOcc) DEALLOCATE (prec_vv) DEALLOCATE (siginvTFTsiginv) DEALLOCATE (STsiginv_0) @@ -1765,6 +1966,17 @@ CONTAINS DEALLOCATE (nocc) DEALLOCATE (m_theta, m_t_in_local) + IF (penalty_occ_local) THEN + DO idim0 = 1, dim_op + DO reim = 1, SIZE(op_sm_set_qs, 1) + DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix) + ENDDO + END DO + DEALLOCATE (op_sm_set_qs) + DEALLOCATE (op_sm_set_almo) + DEALLOCATE (weights) + ENDIF IF (.NOT. converged .AND. .NOT. optimizer%early_stopping_on) THEN CPABORT("Optimization not converged! ") @@ -1774,6 +1986,1006 @@ CONTAINS END SUBROUTINE almo_scf_xalmo_pcg +! ************************************************************************************************** +!> \brief Optimization of NLMOs using PCG minimizers +!> \param qs_env ... +!> \param optimizer controls the optimization algorithm +!> \param matrix_s - AO overlap (NAOs x NAOs) +!> \param matrix_mo_in - initial MOs (NAOs x NMOs) +!> \param matrix_mo_out - final MOs (NAOs x NMOs) +!> \param template_matrix_sigma - template (NMOs x NMOs) +!> \param overlap_determinant - the determinant of the MOs overlap +!> \param mat_distr_aos - info on the distribution of AOs +!> \param virtuals ... +!> \param eps_filter ... +!> \par History +!> 2018.10 created [Rustam Z Khaliullin] +!> \author Rustam Z Khaliullin +! ************************************************************************************************** + SUBROUTINE almo_scf_construct_nlmos(qs_env, optimizer, & + matrix_s, matrix_mo_in, matrix_mo_out, & + template_matrix_sigma, overlap_determinant, & + mat_distr_aos, virtuals, eps_filter) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(optimizer_options_type), INTENT(INOUT) :: optimizer + TYPE(dbcsr_type), INTENT(IN) :: matrix_s + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), & + INTENT(INOUT) :: matrix_mo_in, matrix_mo_out + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), & + INTENT(IN) :: template_matrix_sigma + REAL(KIND=dp), INTENT(INOUT) :: overlap_determinant + INTEGER, INTENT(IN) :: mat_distr_aos + LOGICAL, INTENT(IN) :: virtuals + REAL(KIND=dp), INTENT(IN) :: eps_filter + + CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_construct_nlmos', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=30) :: iter_type, print_string + INTEGER :: cg_iteration, dim_op, handle, iatom, idim0, isgf, ispin, iteration, & + line_search_iteration, linear_search_type, max_iter, natom, ncol, nspins, & + outer_iteration, outer_max_iter, para_group, prec_type, reim, unit_nr + INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nocc, nsgf + LOGICAL :: converged, d_bfgs, just_started, l_bfgs, & + line_search, outer_prepare_to_exit, & + prepare_to_exit, reset_conjugator + REAL(KIND=dp) :: appr_sec_der, beta, bfgs_rho, bfgs_sum, denom, denom2, e0, e1, g0, g0sign, & + g1, g1sign, grad_norm, line_search_error, localization_obj_function, & + localization_obj_function_ispin, next_step_size_guess, obj_function_ispin, objf_diff, & + objf_new, objf_old, penalty_amplitude, penalty_func_ispin, penalty_func_new, spin_factor, & + step_size, t1, t2, tempreal + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diagonal, grad_norm_spin, & + penalty_vol_prefactor, & + suggested_vol_penalty, weights + TYPE(cell_type), POINTER :: cell + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: qs_matrix_s + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set_almo, op_sm_set_qs + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: approx_inv_hessian, bfgs_s, bfgs_y, grad, & + m_S0, m_sig_sqrti_ii, m_siginv, m_sigma, m_t_mo_local, m_theta, m_theta_normalized, & + prev_grad, prev_m_theta, prev_minus_prec_grad, prev_step, step, tempNOcc1, tempOccOcc1, & + tempOccOcc2, tempOccOcc3 + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:, :, :) :: m_B0 + TYPE(lbfgs_history_type) :: nlmo_lbfgs_history + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + + CALL timeset(routineN, handle) + + ! get a useful output_unit + logger => cp_get_default_logger() + IF (logger%para_env%mepos == logger%para_env%source) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + ENDIF + + nspins = SIZE(matrix_mo_in) + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) + IF (.NOT. virtuals) THEN + WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 24), & + " Optimization of occupied NLMOs ", REPEAT("-", 23) + ELSE + WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 24), & + " Optimization of virtual NLMOs ", REPEAT("-", 24) + ENDIF + WRITE (unit_nr, *) + WRITE (unit_nr, '(T2,A13,A6,A23,A14,A14,A9)') "Method", "Iter", & + "Objective Function", "Change", "Convergence", "Time" + WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) + ENDIF + + NULLIFY (particle_set) + + CALL get_qs_env(qs_env=qs_env, & + matrix_s=qs_matrix_s, & + cell=cell, & + particle_set=particle_set, & + qs_kind_set=qs_kind_set) + + natom = SIZE(particle_set, 1) + ALLOCATE (first_sgf(natom)) + ALLOCATE (last_sgf(natom)) + ALLOCATE (nsgf(natom)) + ! construction of + CALL get_particle_set(particle_set, qs_kind_set, & + first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf) + + ! m_theta contains a set of variational parameters + ! that define one-electron orbitals + ALLOCATE (m_theta(nspins)) + DO ispin = 1, nspins + CALL dbcsr_create(m_theta(ispin), & + template=template_matrix_sigma(ispin), & + matrix_type=dbcsr_type_no_symmetry) + ! create initial guess for the main variable - identity matrix + CALL dbcsr_set(m_theta(ispin), 0.0_dp) + CALL dbcsr_add_on_diag(m_theta(ispin), 1.0_dp) + ENDDO + + SELECT CASE (optimizer%opt_penalty%operator_type) + CASE (op_loc_berry) + + IF (cell%orthorhombic) THEN + dim_op = 3 + ELSE + dim_op = 6 + END IF + ALLOCATE (weights(6)) + weights = 0.0_dp + CALL initialize_weights(cell, weights) + ALLOCATE (op_sm_set_qs(2, dim_op)) + ALLOCATE (op_sm_set_almo(2, dim_op)) + ! allocate space for T0^t.B.T0 + ALLOCATE (m_B0(2, dim_op, nspins)) + DO idim0 = 1, dim_op + DO reim = 1, SIZE(op_sm_set_qs, 1) + NULLIFY (op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix) + ALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + ALLOCATE (op_sm_set_almo(reim, idim0)%matrix) + CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, & + name="almo_scf_env%op_sm_"//TRIM(ADJUSTL(cp_to_string(reim)))//"-"//TRIM(ADJUSTL(cp_to_string(idim0)))) + CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp) + CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, matrix_s, & + name="almo_scf_env%op_sm_"//TRIM(ADJUSTL(cp_to_string(reim)))//"-"//TRIM(ADJUSTL(cp_to_string(idim0)))) + CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp) + DO ispin = 1, nspins + CALL dbcsr_create(m_B0(reim, idim0, ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_set(m_B0(reim, idim0, ispin), 0.0_dp) + ENDDO + ENDDO + ENDDO + + CALL compute_berry_operator(qs_env, cell, op_sm_set_qs, dim_op) + + CASE (op_loc_pipek) + + dim_op = natom + ALLOCATE (weights(dim_op)) + weights = 1.0_dp + + ALLOCATE (m_B0(1, dim_op, nspins)) + !m_B0 first dim is 1 now! + DO idim0 = 1, dim_op + DO reim = 1, 1 !SIZE(op_sm_set_qs, 1) + DO ispin = 1, nspins + CALL dbcsr_create(m_B0(reim, idim0, ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_set(m_B0(reim, idim0, ispin), 0.0_dp) + ENDDO + ENDDO + ENDDO + + END SELECT + + ! penalty amplitude adjusts the strenght of volume conservation + penalty_amplitude = optimizer%opt_penalty%penalty_strength + !penalty_occ_vol = ( optimizer%opt_penalty%occ_vol_method .NE. penalty_type_none ) + !penalty_local = ( optimizer%opt_penalty%occ_loc_method .NE. penalty_type_none ) + + ! preconditioner control + prec_type = optimizer%preconditioner + + ! use diagonal BFGS if preconditioner is set + d_bfgs = .FALSE. + l_bfgs = .FALSE. + IF (prec_type .NE. xalmo_prec_zero) l_bfgs = .TRUE. + IF (l_bfgs .AND. (optimizer%conjugator .NE. cg_zero)) THEN + CPABORT("Cannot use conjugators with BFGS") + ENDIF + IF (l_bfgs) THEN + CALL lbfgs_create(nlmo_lbfgs_history, nspins, nstore=10) + ENDIF + + IF (nspins == 1) THEN + spin_factor = 2.0_dp + ELSE + spin_factor = 1.0_dp + ENDIF + + ALLOCATE (grad_norm_spin(nspins)) + ALLOCATE (nocc(nspins)) + ALLOCATE (penalty_vol_prefactor(nspins)) + ALLOCATE (suggested_vol_penalty(nspins)) + + ! create a local copy of matrix_mo_in because + ! matrix_mo_in and matrix_mo_out can be the same matrix + ! we need to make sure data in matrix_mo_in is intact + ! after we start writing to matrix_mo_out + ALLOCATE (m_t_mo_local(nspins)) + DO ispin = 1, nspins + CALL dbcsr_create(m_t_mo_local(ispin), & + template=matrix_mo_in(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_copy(m_t_mo_local(ispin), matrix_mo_in(ispin)) + ENDDO + + ALLOCATE (approx_inv_hessian(nspins)) + ALLOCATE (m_theta_normalized(nspins)) + ALLOCATE (prev_m_theta(nspins)) + ALLOCATE (m_S0(nspins)) + ALLOCATE (prev_grad(nspins)) + ALLOCATE (grad(nspins)) + ALLOCATE (prev_step(nspins)) + ALLOCATE (step(nspins)) + ALLOCATE (prev_minus_prec_grad(nspins)) + ALLOCATE (m_sig_sqrti_ii(nspins)) + ALLOCATE (m_sigma(nspins)) + ALLOCATE (m_siginv(nspins)) + ALLOCATE (tempNOcc1(nspins)) + ALLOCATE (tempOccOcc1(nspins)) + ALLOCATE (tempOccOcc2(nspins)) + ALLOCATE (tempOccOcc3(nspins)) + ALLOCATE (bfgs_y(nspins)) + ALLOCATE (bfgs_s(nspins)) + + DO ispin = 1, nspins + + ! init temporary storage + CALL dbcsr_create(tempNOcc1(ispin), & + template=matrix_mo_out(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(approx_inv_hessian(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_theta_normalized(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(prev_m_theta(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_S0(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(prev_grad(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(grad(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(prev_step(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(step(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(prev_minus_prec_grad(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_sig_sqrti_ii(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_sigma(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_siginv(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc1(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc2(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc3(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(bfgs_s(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(bfgs_y(ispin), & + template=m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + + CALL dbcsr_set(step(ispin), 0.0_dp) + CALL dbcsr_set(prev_step(ispin), 0.0_dp) + + CALL dbcsr_get_info(template_matrix_sigma(ispin), & + nfullrows_total=nocc(ispin)) + + penalty_vol_prefactor(ispin) = -penalty_amplitude !KEEP: * spin_factor * nocc(ispin) + + ! compute m_S0=T0^t.S.T0 + CALL dbcsr_multiply("N", "N", 1.0_dp, & + matrix_s, & + m_t_mo_local(ispin), & + 0.0_dp, tempNOcc1(ispin), & + filter_eps=eps_filter) + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_t_mo_local(ispin), & + tempNOcc1(ispin), & + 0.0_dp, m_S0(ispin), & + filter_eps=eps_filter) + + SELECT CASE (optimizer%opt_penalty%operator_type) + + CASE (op_loc_berry) + + ! compute m_B0=T0^t.B.T0 + DO idim0 = 1, SIZE(op_sm_set_qs, 2) ! this loop is over miller ind + + DO reim = 1, SIZE(op_sm_set_qs, 1) ! this loop is over Re/Im + + CALL matrix_qs_to_almo(op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix, & + mat_distr_aos, .FALSE.) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + op_sm_set_almo(reim, idim0)%matrix, & + m_t_mo_local(ispin), & + 0.0_dp, tempNOcc1(ispin), & + filter_eps=eps_filter) + + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_t_mo_local(ispin), & + tempNOcc1(ispin), & + 0.0_dp, m_B0(reim, idim0, ispin), & + filter_eps=eps_filter) + + DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix) + + ENDDO + + ENDDO ! end loop over idim0 + + CASE (op_loc_pipek) + + ! compute m_B0=T0^t.B.T0 + DO iatom = 1, natom ! this loop is over "miller" ind + + isgf = first_sgf(iatom) + ncol = nsgf(iatom) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + matrix_s, & + m_t_mo_local(ispin), & + 0.0_dp, tempNOcc1(ispin), & + filter_eps=eps_filter) + + CALL dbcsr_multiply("T", "N", 0.5_dp, & + m_t_mo_local(ispin), & + tempNOcc1(ispin), & + 0.0_dp, m_B0(1, iatom, ispin), & + first_k=isgf, last_k=isgf + ncol - 1, & + filter_eps=eps_filter) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + matrix_s, & + m_t_mo_local(ispin), & + 0.0_dp, tempNOcc1(ispin), & + first_k=isgf, last_k=isgf + ncol - 1, & + filter_eps=eps_filter) + + CALL dbcsr_multiply("T", "N", 0.5_dp, & + m_t_mo_local(ispin), & + tempNOcc1(ispin), & + 1.0_dp, m_B0(1, iatom, ispin), & + filter_eps=eps_filter) + + ENDDO ! end loop over iatom + + END SELECT + + ENDDO ! ispin + + IF (optimizer%opt_penalty%operator_type .EQ. op_loc_berry) THEN + DO idim0 = 1, SIZE(op_sm_set_qs, 2) ! this loop is over miller ind + DO reim = 1, SIZE(op_sm_set_qs, 1) ! this loop is over Re/Im + DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix) + ENDDO + ENDDO + DEALLOCATE (op_sm_set_qs, op_sm_set_almo) + ENDIF + + ! start the outer SCF loop + outer_max_iter = optimizer%max_iter_outer_loop + outer_prepare_to_exit = .FALSE. + outer_iteration = 0 + grad_norm = 0.0_dp + penalty_func_new = 0.0_dp + linear_search_type = 1 ! safe restart, no quadratic assumption, takes more steps + localization_obj_function = 0.0_dp + penalty_func_new = 0.0_dp + + DO + + ! start the inner SCF loop + max_iter = optimizer%max_iter + prepare_to_exit = .FALSE. + line_search = .FALSE. + converged = .FALSE. + iteration = 0 + cg_iteration = 0 + line_search_iteration = 0 + obj_function_ispin = 0.0_dp + objf_new = 0.0_dp + objf_old = 0.0_dp + objf_diff = 0.0_dp + line_search_error = 0.0_dp + t1 = m_walltime() + next_step_size_guess = 0.0_dp + + DO + + just_started = (iteration .EQ. 0) .AND. (outer_iteration .EQ. 0) + + DO ispin = 1, nspins + + CALL dbcsr_get_info(m_sig_sqrti_ii(ispin), group=para_group) + + ! compute diagonal (a^t.sigma0.a)^(-1/2) + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_S0(ispin), m_theta(ispin), 0.0_dp, & + tempOccOcc1(ispin), & + filter_eps=eps_filter) + CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp) + CALL dbcsr_add_on_diag(m_sig_sqrti_ii(ispin), 1.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_theta(ispin), tempOccOcc1(ispin), 0.0_dp, & + m_sig_sqrti_ii(ispin), & + retain_sparsity=.TRUE.) + ALLOCATE (diagonal(nocc(ispin))) + CALL dbcsr_get_diag(m_sig_sqrti_ii(ispin), diagonal) + CALL mp_sum(diagonal, para_group) + ! TODO: works for zero diagonal elements? + diagonal(:) = 1.0_dp/SQRT(diagonal(:)) + CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp) + CALL dbcsr_set_diag(m_sig_sqrti_ii(ispin), diagonal) + DEALLOCATE (diagonal) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_theta(ispin), & + m_sig_sqrti_ii(ispin), & + 0.0_dp, m_theta_normalized(ispin), & + filter_eps=eps_filter) + + ! compute new orbitals + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_t_mo_local(ispin), & + m_theta_normalized(ispin), & + 0.0_dp, matrix_mo_out(ispin), & + filter_eps=eps_filter) + + ENDDO + + ! compute objective function + localization_obj_function = 0.0_dp + penalty_func_new = 0.0_dp + DO ispin = 1, nspins + + CALL compute_obj_nlmos( & + !obj_function_ispin=obj_function_ispin, & + localization_obj_function_ispin=localization_obj_function_ispin, & + penalty_func_ispin=penalty_func_ispin, & + overlap_determinant=overlap_determinant, & + m_sigma=m_sigma(ispin), & + nocc=nocc(ispin), & + m_B0=m_B0(:, :, ispin), & + m_theta_normalized=m_theta_normalized(ispin), & + template_matrix_mo=matrix_mo_out(ispin), & + weights=weights, & + m_S0=m_S0(ispin), & + just_started=just_started, & + penalty_vol_prefactor=penalty_vol_prefactor(ispin), & + penalty_amplitude=penalty_amplitude, & + eps_filter=eps_filter) + + localization_obj_function = localization_obj_function + localization_obj_function_ispin + penalty_func_new = penalty_func_new + penalty_func_ispin + + ENDDO ! ispin + objf_new = penalty_func_new + localization_obj_function + + DO ispin = 1, nspins + ! save the previous gradient to compute beta + ! do it only if the previous grad was computed + ! for .NOT.line_search + IF (line_search_iteration .EQ. 0 .AND. iteration .NE. 0) THEN + CALL dbcsr_copy(prev_grad(ispin), grad(ispin)) + ENDIF + + ENDDO ! ispin + + ! compute the gradient + DO ispin = 1, nspins + + CALL invert_Hotelling( & + matrix_inverse=m_siginv(ispin), & + matrix=m_sigma(ispin), & + threshold=eps_filter*10.0_dp, & + filter_eps=eps_filter, & + silent=.FALSE.) + + CALL compute_gradient_nlmos( & + m_grad_out=grad(ispin), & + m_B0=m_B0(:, :, ispin), & + weights=weights, & + m_S0=m_S0(ispin), & + m_theta_normalized=m_theta_normalized(ispin), & + m_siginv=m_siginv(ispin), & + m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), & + penalty_vol_prefactor=penalty_vol_prefactor(ispin), & + eps_filter=eps_filter, & + suggested_vol_penalty=suggested_vol_penalty(ispin)) + + ENDDO ! ispin + + ! check convergence and other exit criteria + DO ispin = 1, nspins + CALL dbcsr_norm(grad(ispin), dbcsr_norm_maxabsnorm, & + norm_scalar=grad_norm_spin(ispin)) + ENDDO ! ispin + grad_norm = MAXVAL(grad_norm_spin) + + converged = (grad_norm .LE. optimizer%eps_error) + IF (converged .OR. (iteration .GE. max_iter)) THEN + prepare_to_exit = .TRUE. + ENDIF + + ! it is not time to exit just yet + IF (.NOT. prepare_to_exit) THEN + + ! check the gradient along the step direction + ! and decide whether to switch to the line-search mode + ! do not do this in the first iteration + IF (iteration .NE. 0) THEN + + ! enforce at least one line search + ! without even checking the error + IF (.NOT. line_search) THEN + + line_search = .TRUE. + line_search_iteration = line_search_iteration + 1 + + ELSE + + ! check the line-search error and decide whether to + ! change the direction + line_search_error = 0.0_dp + denom = 0.0_dp + denom2 = 0.0_dp + + DO ispin = 1, nspins + + CALL dbcsr_dot(grad(ispin), step(ispin), tempreal) + line_search_error = line_search_error + tempreal + CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal) + denom = denom + tempreal + CALL dbcsr_dot(step(ispin), step(ispin), tempreal) + denom2 = denom2 + tempreal + + ENDDO ! ispin + + ! cosine of the angle between the step and grad + ! (must be close to zero at convergence) + line_search_error = line_search_error/SQRT(denom)/SQRT(denom2) + + IF (ABS(line_search_error) .GT. optimizer%lin_search_eps_error) THEN + line_search = .TRUE. + line_search_iteration = line_search_iteration + 1 + ELSE + line_search = .FALSE. + line_search_iteration = 0 + ENDIF + + ENDIF + + ENDIF ! iteration.ne.0 + + IF (line_search) THEN + objf_diff = 0.0_dp + ELSE + objf_diff = objf_new - objf_old + objf_old = objf_new + ENDIF + + ! update the step direction + IF (.NOT. line_search) THEN + + cg_iteration = cg_iteration + 1 + + ! save the previous step + DO ispin = 1, nspins + CALL dbcsr_copy(prev_step(ispin), step(ispin)) + ENDDO ! ispin + + ! compute the new step: + ! if available use second derivative info - bfgs, hessian, preconditioner + IF (prec_type .EQ. xalmo_prec_zero) THEN ! no second derivatives + + ! no preconditioner + DO ispin = 1, nspins + + CALL dbcsr_copy(step(ispin), grad(ispin)) + CALL dbcsr_scale(step(ispin), -1.0_dp) + + ENDDO ! ispin + + ELSE ! use second derivatives + + ! compute and invert hessian/precond? + IF (iteration .EQ. 0) THEN + + IF (d_bfgs) THEN + + ! create matrix filled with 1.0 here + CALL fill_matrix_with_ones(approx_inv_hessian(1)) + IF (nspins .GT. 1) THEN + DO ispin = 2, nspins + CALL dbcsr_copy(approx_inv_hessian(ispin), approx_inv_hessian(1)) + ENDDO + ENDIF + + ELSE IF (l_bfgs) THEN + + CALL lbfgs_seed(nlmo_lbfgs_history, m_theta, grad) + DO ispin = 1, nspins + CALL dbcsr_copy(step(ispin), grad(ispin)) + CALL dbcsr_scale(step(ispin), -1.0_dp) + ENDDO ! ispin + + ELSE + + ! computing preconditioner + DO ispin = 1, nspins + + ! TODO: write preconditioner code later + ! For now, create matrix filled with 1.0 here + CALL fill_matrix_with_ones(approx_inv_hessian(ispin)) + !CALL compute_preconditioner(& + ! m_prec_out=approx_hessian(ispin),& + ! m_ks=almo_scf_env%matrix_ks(ispin),& + ! m_s=matrix_s,& + ! m_siginv=almo_scf_env%template_matrix_sigma(ispin),& + ! m_quench_t=quench_t(ispin),& + ! m_FTsiginv=FTsiginv(ispin),& + ! m_siginvTFTsiginv=siginvTFTsiginv(ispin),& + ! m_ST=ST(ispin),& + ! para_env=almo_scf_env%para_env,& + ! blacs_env=almo_scf_env%blacs_env,& + ! nocc_of_domain=almo_scf_env%nocc_of_domain(:,ispin),& + ! domain_s_inv=almo_scf_env%domain_s_inv(:,ispin),& + ! domain_r_down=domain_r_down(:,ispin),& + ! cpu_of_domain=almo_scf_env%cpu_of_domain,& + ! domain_map=almo_scf_env%domain_map(ispin),& + ! assume_t0_q0x=assume_t0_q0x,& + ! penalty_occ_vol=penalty_occ_vol,& + ! penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin),& + ! eps_filter=eps_filter,& + ! neg_thr=0.5_dp,& + ! spin_factor=spin_factor,& + ! special_case=my_special_case) + !CALL invert hessian + ENDDO ! ispin + + ENDIF + + ELSE ! not iteration zero + + ! update approx inverse hessian + IF (d_bfgs) THEN ! diagonal BFGS + + DO ispin = 1, nspins + + ! compute s and y + CALL dbcsr_copy(bfgs_y(ispin), grad(ispin)) + CALL dbcsr_add(bfgs_y(ispin), prev_grad(ispin), 1.0_dp, -1.0_dp) + CALL dbcsr_copy(bfgs_s(ispin), m_theta(ispin)) + CALL dbcsr_add(bfgs_s(ispin), prev_m_theta(ispin), 1.0_dp, -1.0_dp) + + ! compute rho + CALL dbcsr_dot(grad(ispin), step(ispin), bfgs_rho) + bfgs_rho = 1.0_dp/bfgs_rho + + ! compute the sum of the squared elements of bfgs_y + CALL dbcsr_dot(bfgs_y(ispin), bfgs_y(ispin), bfgs_sum) + + ! first term: start collecting new inv hessian in this temp matrix + CALL dbcsr_copy(tempOccOcc2(ispin), approx_inv_hessian(ispin)) + + ! second term: + rho * s * s + CALL dbcsr_hadamard_product(bfgs_s(ispin), bfgs_s(ispin), tempOccOcc1(ispin)) + CALL dbcsr_add(tempOccOcc2(ispin), tempOccOcc1(ispin), 1.0_dp, bfgs_rho) + + ! third term: + rho^2 * s * s * H * sum_(y * y) + CALL dbcsr_hadamard_product(tempOccOcc1(ispin), & + approx_inv_hessian(ispin), tempOccOcc3(ispin)) + CALL dbcsr_add(tempOccOcc2(ispin), tempOccOcc3(ispin), & + 1.0_dp, bfgs_rho*bfgs_rho*bfgs_sum) + + ! fourth term: - 2 * rho * s * y * H + CALL dbcsr_hadamard_product(bfgs_y(ispin), & + approx_inv_hessian(ispin), tempOccOcc1(ispin)) + CALL dbcsr_hadamard_product(bfgs_s(ispin), tempOccOcc1(ispin), tempOccOcc3(ispin)) + CALL dbcsr_add(tempOccOcc2(ispin), tempOccOcc3(ispin), & + 1.0_dp, -2.0_dp*bfgs_rho) + + CALL dbcsr_copy(approx_inv_hessian(ispin), tempOccOcc2(ispin)) + + ENDDO + + ELSE IF (l_bfgs) THEN + + CALL lbfgs_get_direction(nlmo_lbfgs_history, m_theta, grad, step) + + ENDIF ! which method? + + ENDIF ! compute approximate inverse hessian + + IF (.NOT. l_bfgs) THEN + + DO ispin = 1, nspins + + CALL dbcsr_hadamard_product(approx_inv_hessian(ispin), & + grad(ispin), step(ispin)) + CALL dbcsr_scale(step(ispin), -1.0_dp) + + ENDDO ! ispin + + ENDIF + + ENDIF ! second derivative type fork + + ! check whether we need to reset conjugate directions + IF (iteration .EQ. 0) THEN + reset_conjugator = .TRUE. + ENDIF + + ! compute the conjugation coefficient - beta + IF (.NOT. reset_conjugator) THEN + CALL compute_cg_beta( & + beta=beta, & + reset_conjugator=reset_conjugator, & + conjugator=optimizer%conjugator, & + grad=grad(:), & + prev_grad=prev_grad(:), & + step=step(:), & + prev_step=prev_step(:), & + prev_minus_prec_grad=prev_minus_prec_grad(:) & + ) + + ENDIF + + IF (reset_conjugator) THEN + + beta = 0.0_dp + IF (unit_nr > 0 .AND. (.NOT. just_started)) THEN + WRITE (unit_nr, '(T2,A35)') "Re-setting conjugator to zero" + ENDIF + reset_conjugator = .FALSE. + + ENDIF + + ! save the preconditioned gradient (useful for beta) + DO ispin = 1, nspins + + CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin)) + + ! conjugate the step direction + CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta) + + ENDDO ! ispin + + ENDIF ! update the step direction + + ! estimate the step size + IF (.NOT. line_search) THEN + ! we just changed the direction and + ! we have only E and grad from the current step + ! it is not enough to compute step_size - just guess it + e0 = objf_new + g0 = 0.0_dp + DO ispin = 1, nspins + CALL dbcsr_dot(grad(ispin), step(ispin), tempreal) + g0 = g0 + tempreal + ENDDO ! ispin + g0sign = SIGN(1.0_dp, g0) ! sign of g0 + IF (linear_search_type .EQ. 1) THEN ! this is quadratic LS + IF (iteration .EQ. 0) THEN + step_size = optimizer%lin_search_step_size_guess + ELSE + IF (next_step_size_guess .LE. 0.0_dp) THEN + step_size = optimizer%lin_search_step_size_guess + ELSE + ! take the last value + step_size = optimizer%lin_search_step_size_guess + !step_size = next_step_size_guess*1.05_dp + ENDIF + ENDIF + ELSE IF (linear_search_type .EQ. 2) THEN ! this is cautious LS + ! this LS type is designed not to trust quadratic appr + ! so it always restarts from a safe step size + step_size = optimizer%lin_search_step_size_guess + ENDIF + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T21,3A19)') "Line position", "Line grad", "Next line step" + WRITE (unit_nr, '(T2,A19,3F19.5)') "Line search", 0.0_dp, g0, step_size + ENDIF + next_step_size_guess = step_size + ELSE ! this is not the first line search + e1 = objf_new + g1 = 0.0_dp + DO ispin = 1, nspins + CALL dbcsr_dot(grad(ispin), step(ispin), tempreal) + g1 = g1 + tempreal + ENDDO ! ispin + g1sign = SIGN(1.0_dp, g1) ! sign of g1 + IF (linear_search_type .EQ. 1) THEN + ! we have accumulated some points along this direction + ! use only the most recent g0 (quadratic approximation) + appr_sec_der = (g1 - g0)/step_size + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, '(A2,7F12.5)') & + ! "DT", e0, e1, g0, g1, appr_sec_der, step_size, -g1/appr_sec_der + !ENDIF + step_size = -g1/appr_sec_der + ELSE IF (linear_search_type .EQ. 2) THEN + ! alternative method for finding step size + ! do not use quadratic approximation, only gradient signs + IF (g1sign .NE. g0sign) THEN + step_size = -step_size/2.0; + ELSE + step_size = step_size*1.5; + ENDIF + ENDIF + ! end alternative LS types + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T21,3A19)') "Line position", "Line grad", "Next line step" + WRITE (unit_nr, '(T2,A19,3F19.5)') "Line search", next_step_size_guess, g1, step_size + ENDIF + e0 = e1 + g0 = g1 + g0sign = g1sign + next_step_size_guess = next_step_size_guess + step_size + ENDIF + + ! update theta + DO ispin = 1, nspins + IF (.NOT. line_search) THEN ! we prepared to perform the first line search + ! "previous" refers to the previous CG step, not the previous LS step + CALL dbcsr_copy(prev_m_theta(ispin), m_theta(ispin)) + ENDIF + CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size) + ENDDO ! ispin + + ENDIF ! not.prepare_to_exit + + IF (line_search) THEN + iter_type = "LS" + ELSE + iter_type = "CG" + ENDIF + + t2 = m_walltime() + IF (unit_nr > 0) THEN + iter_type = TRIM("NLMO OPT "//iter_type) + WRITE (unit_nr, '(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') & + iter_type, iteration, & + objf_new, objf_diff, grad_norm, & + t2 - t1 + WRITE (unit_nr, '(T2,A19,F23.10)') & + "Localization:", localization_obj_function + WRITE (unit_nr, '(T2,A19,F23.10)') & + "Orthogonalization:", penalty_func_new + ENDIF + t1 = m_walltime() + + iteration = iteration + 1 + IF (prepare_to_exit) EXIT + + ENDDO ! inner loop + + IF (converged .OR. (outer_iteration .GE. outer_max_iter)) THEN + outer_prepare_to_exit = .TRUE. + ENDIF + + outer_iteration = outer_iteration + 1 + IF (outer_prepare_to_exit) EXIT + + ENDDO ! outer loop + + ! return the optimal determinant penalty + optimizer%opt_penalty%penalty_strength = 0.0_dp + DO ispin = 1, nspins + optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength + & + (-1.0_dp)*penalty_vol_prefactor(ispin) + ENDDO + optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength/nspins + + IF (converged) THEN + iter_type = "Final" + ELSE + iter_type = "Unconverged" + ENDIF + + IF (unit_nr > 0) THEN + WRITE (unit_nr, '()') + print_string = TRIM(iter_type)//" localization:" + WRITE (unit_nr, '(T2,A29,F30.10)') & + print_string, localization_obj_function + print_string = TRIM(iter_type)//" determinant:" + WRITE (unit_nr, '(T2,A29,F30.10)') & + print_string, overlap_determinant + print_string = TRIM(iter_type)//" penalty strength:" + WRITE (unit_nr, '(T2,A29,F30.10)') & + print_string, optimizer%opt_penalty%penalty_strength + ENDIF + + ! clean up + IF (l_bfgs) THEN + CALL lbfgs_release(nlmo_lbfgs_history) + ENDIF + DO ispin = 1, nspins + DO idim0 = 1, SIZE(m_B0, 2) + DO reim = 1, SIZE(m_B0, 1) + CALL dbcsr_release(m_B0(reim, idim0, ispin)) + ENDDO + ENDDO + CALL dbcsr_release(m_theta(ispin)) + CALL dbcsr_release(m_t_mo_local(ispin)) + CALL dbcsr_release(tempNOcc1(ispin)) + CALL dbcsr_release(approx_inv_hessian(ispin)) + CALL dbcsr_release(prev_m_theta(ispin)) + CALL dbcsr_release(m_theta_normalized(ispin)) + CALL dbcsr_release(m_S0(ispin)) + CALL dbcsr_release(prev_grad(ispin)) + CALL dbcsr_release(grad(ispin)) + CALL dbcsr_release(prev_step(ispin)) + CALL dbcsr_release(step(ispin)) + CALL dbcsr_release(prev_minus_prec_grad(ispin)) + CALL dbcsr_release(m_sig_sqrti_ii(ispin)) + CALL dbcsr_release(m_sigma(ispin)) + CALL dbcsr_release(m_siginv(ispin)) + CALL dbcsr_release(tempOccOcc1(ispin)) + CALL dbcsr_release(tempOccOcc2(ispin)) + CALL dbcsr_release(tempOccOcc3(ispin)) + CALL dbcsr_release(bfgs_y(ispin)) + CALL dbcsr_release(bfgs_s(ispin)) + ENDDO ! ispin + + DEALLOCATE (grad_norm_spin) + DEALLOCATE (nocc) + DEALLOCATE (penalty_vol_prefactor) + DEALLOCATE (suggested_vol_penalty) + + DEALLOCATE (approx_inv_hessian) + DEALLOCATE (prev_m_theta) + DEALLOCATE (m_theta_normalized) + DEALLOCATE (m_S0) + DEALLOCATE (prev_grad) + DEALLOCATE (grad) + DEALLOCATE (prev_step) + DEALLOCATE (step) + DEALLOCATE (prev_minus_prec_grad) + DEALLOCATE (m_sig_sqrti_ii) + DEALLOCATE (m_sigma) + DEALLOCATE (m_siginv) + DEALLOCATE (tempNOcc1) + DEALLOCATE (tempOccOcc1) + DEALLOCATE (tempOccOcc2) + DEALLOCATE (tempOccOcc3) + DEALLOCATE (bfgs_y) + DEALLOCATE (bfgs_s) + + DEALLOCATE (m_theta, m_t_mo_local) + DEALLOCATE (m_B0) + DEALLOCATE (weights) + DEALLOCATE (first_sgf, last_sgf, nsgf) + + IF (.NOT. converged) THEN + CPABORT("Optimization not converged! ") + ENDIF + + CALL timestop(handle) + + END SUBROUTINE almo_scf_construct_nlmos + ! ************************************************************************************************** !> \brief Analysis of the orbitals !> \param detailed_analysis ... @@ -2010,7 +3222,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE split_v_blk(almo_scf_env) - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'split_v_blk', routineP = moduleN//':'//routineN @@ -2079,7 +3291,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE harris_foulkes_correction(almo_scf_env) - TYPE(almo_scf_env_type) :: almo_scf_env + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env CHARACTER(len=*), PARAMETER :: routineN = 'harris_foulkes_correction', & routineP = moduleN//':'//routineN @@ -2112,18 +3324,6 @@ CONTAINS vr_index_sqrt_inv TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_p_almo_scf_converged -!TYPE(dbcsr_iterator_type) :: iter -!TYPE(dbcsr_type) :: tmp11,tmp22,tmp33 -!REAL(kind=dp) :: k_var1, k_var2 -!TYPE(dbcsr_type) :: fake_step -! -!TYPE(dbcsr_type) :: sigma_dr, sigma_dr2, sigma_rr, sigma_rr2 -!TYPE(dbcsr_type) :: fake_a,fake_b,fake_k0 -!INTEGER :: retained_v,discarded_v,i_row,j_col -!TYPE(dbcsr_type) :: matrix_rst0, matrix_rst1, matrix_rst2, ss_vv -!REAL(KIND=dp) :: filter_memorize, init_filter, occ_vv -!INTEGER :: ppp - CALL timeset(routineN, handle) ! get a useful output_unit @@ -4737,12 +5937,18 @@ CONTAINS !> \param optimize_theta ... !> \param normalize_orbitals ... !> \param penalty_occ_vol ... +!> \param penalty_occ_local ... !> \param penalty_occ_vol_prefactor ... !> \param envelope_amplitude ... !> \param eps_filter ... !> \param spin_factor ... !> \param special_case ... !> \param m_sig_sqrti_ii ... +!> \param op_sm_set ... +!> \param weights ... +!> \param energy_coeff ... +!> \param localiz_coeff ... +!> \param z2 ... !> \par History !> 2015.03 created [Rustam Z Khaliullin] !> \author Rustam Z Khaliullin @@ -4751,9 +5957,10 @@ CONTAINS m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv0, & m_theta, domain_s_inv, domain_r_down, & cpu_of_domain, domain_map, assume_t0_q0x, optimize_theta, & - normalize_orbitals, penalty_occ_vol, & + normalize_orbitals, penalty_occ_vol, penalty_occ_local, & penalty_occ_vol_prefactor, envelope_amplitude, eps_filter, spin_factor, & - special_case, m_sig_sqrti_ii) + special_case, m_sig_sqrti_ii, op_sm_set, weights, energy_coeff, & + localiz_coeff, z2) TYPE(dbcsr_type), INTENT(INOUT) :: m_grad_out TYPE(dbcsr_type), INTENT(IN) :: m_ks, m_s, m_t, m_t0, m_siginv, & @@ -4766,21 +5973,30 @@ CONTAINS TYPE(domain_map_type), INTENT(IN) :: domain_map LOGICAL, INTENT(IN) :: assume_t0_q0x, optimize_theta, & normalize_orbitals, penalty_occ_vol + LOGICAL, INTENT(IN), OPTIONAL :: penalty_occ_local REAL(KIND=dp), INTENT(IN) :: penalty_occ_vol_prefactor, & envelope_amplitude, eps_filter, & spin_factor INTEGER, INTENT(IN) :: special_case TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: m_sig_sqrti_ii + TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & + POINTER :: op_sm_set + REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: weights + REAL(KIND=dp), INTENT(IN), OPTIONAL :: energy_coeff, localiz_coeff + REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: z2 CHARACTER(len=*), PARAMETER :: routineN = 'compute_gradient', & routineP = moduleN//':'//routineN - INTEGER :: handle, nao - REAL(KIND=dp) :: energy_g_norm, penalty_occ_vol_g_norm + INTEGER :: dim0, handle, idim0, ielem, nao, reim + LOGICAL :: my_penalty_local + REAL(KIND=dp) :: coeff, energy_g_norm, my_energy_coeff, & + my_localiz_coeff, & + penalty_occ_vol_g_norm + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tg_diagonal TYPE(dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_no_3, & - m_tmp_oo_1, m_tmp_oo_2 - -!TYPE(dbcsr_type), INTENT(INOUT) :: m_tmp_no_2 + m_tmp_oo_1, m_tmp_oo_2, temp1, temp2, & + tempNOcc1, tempOccOcc1 CALL timeset(routineN, handle) @@ -4788,6 +6004,19 @@ CONTAINS CPABORT("Normalization matrix is required") ENDIF + my_penalty_local = .FALSE. + my_localiz_coeff = 1.0_dp + my_energy_coeff = 0.0_dp + IF (PRESENT(localiz_coeff)) THEN + my_localiz_coeff = localiz_coeff + ENDIF + IF (PRESENT(energy_coeff)) THEN + my_energy_coeff = energy_coeff + ENDIF + IF (PRESENT(penalty_occ_local)) THEN + my_penalty_local = penalty_occ_local + ENDIF + ! use this otherways unused variables CALL dbcsr_get_info(matrix=m_ks, nfullrows_total=nao) CALL dbcsr_get_info(matrix=m_s, nfullrows_total=nao) @@ -4808,6 +6037,18 @@ CONTAINS CALL dbcsr_create(m_tmp_oo_2, & template=m_siginv, & matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempNOcc1, & + template=m_t, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc1, & + template=m_siginv, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(temp1, & + template=m_t, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(temp2, & + template=m_t, & + matrix_type=dbcsr_type_no_symmetry) ! do d_E/d_T first !IF (.NOT.PRESENT(m_FTsiginv)) THEN @@ -4854,6 +6095,58 @@ CONTAINS retain_sparsity=.TRUE.) CALL dbcsr_scale(m_tmp_no_2, 2.0_dp*spin_factor) + ! LzL Add gradient for Localization + IF (my_penalty_local) THEN + + CALL dbcsr_set(temp2, 0.0_dp) ! accumulate the localization gradient here + + DO idim0 = 1, SIZE(op_sm_set, 2) ! this loop is over miller ind + + DO reim = 1, SIZE(op_sm_set, 1) ! this loop is over Re/Im + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + op_sm_set(reim, idim0)%matrix, & + m_t, & + 0.0_dp, tempNOcc1, & + filter_eps=eps_filter) + + ! warning - save time by computing only the diagonal elements + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_t, & + tempNOcc1, & + 0.0_dp, tempOccOcc1, & + filter_eps=eps_filter) + + CALL dbcsr_get_info(tempOccOcc1, nfullrows_total=dim0) + ALLOCATE (tg_diagonal(dim0)) + CALL dbcsr_get_diag(tempOccOcc1, tg_diagonal) + CALL dbcsr_set(tempOccOcc1, 0.0_dp) + CALL dbcsr_set_diag(tempOccOcc1, tg_diagonal) + DEALLOCATE (tg_diagonal) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + tempNOcc1, & + tempOccOcc1, & + 0.0_dp, temp1, & + filter_eps=eps_filter) + + ENDDO + + SELECT CASE (2) ! allows for selection of different spread functionals + CASE (1) ! functional = -W_I * log( |z_I|^2 ) + coeff = -(weights(idim0)/z2(ielem)) + CASE (2) ! functional = W_I * ( 1 - |z_I|^2 ) + coeff = -weights(idim0) + CASE (3) ! functional = W_I * ( 1 - |z_I| ) + coeff = -(weights(idim0)/(2.0_dp*z2(ielem))) + END SELECT + CALL dbcsr_add(temp2, temp1, 1.0_dp, coeff) + !CALL dbcsr_add(grad_loc, temp1, 1.0_dp, 1.0_dp) + + ENDDO ! end loop over idim0 + CALL dbcsr_add(m_tmp_no_2, temp2, my_energy_coeff, my_localiz_coeff*4.0_dp) + ENDIF + ! add penalty on the occupied volume: det(sigma) IF (penalty_occ_vol) THEN !RZK-warning CALL dbcsr_multiply("N","N",& @@ -4870,7 +6163,7 @@ CONTAINS m_siginv, & 0.0_dp, m_tmp_no_1, & retain_sparsity=.TRUE.) -! this norm does not contain the normalization factors + ! this norm does not contain the normalization factors CALL dbcsr_norm(m_tmp_no_1, dbcsr_norm_maxabsnorm, & norm_scalar=penalty_occ_vol_g_norm) CALL dbcsr_norm(m_tmp_no_2, dbcsr_norm_maxabsnorm, & @@ -4888,9 +6181,19 @@ CONTAINS ! where c0 = penalty_occ_vol_prefactor ! This is because tr(T).G_Energy = 0 and ! tr(T).G_Penalty = c0*I - CALL dbcsr_copy(m_tmp_no_1, m_quench_t) - CALL dbcsr_copy(m_tmp_no_1, m_ST, keep_sparsity=.TRUE.) - CALL dbcsr_add(m_tmp_no_2, m_tmp_no_1, 1.0_dp, -penalty_occ_vol_prefactor) + + !! faster way to take the norm into account (tested for vol penalty olny) + !!CALL dbcsr_copy(m_tmp_no_1, m_quench_t) + !!CALL dbcsr_copy(m_tmp_no_1, m_ST, keep_sparsity=.TRUE.) + !!CALL dbcsr_add(m_tmp_no_2, m_tmp_no_1, 1.0_dp, -penalty_occ_vol_prefactor) + !!CALL dbcsr_copy(m_tmp_no_1, m_quench_t) + !!CALL dbcsr_multiply("N", "N", 1.0_dp, & + !! m_tmp_no_2, & + !! m_sig_sqrti_ii, & + !! 0.0_dp, m_tmp_no_1, & + !! retain_sparsity=.TRUE.) + + ! slower way of taking the norm into account CALL dbcsr_copy(m_tmp_no_1, m_quench_t) CALL dbcsr_multiply("N", "N", 1.0_dp, & m_tmp_no_2, & @@ -4898,42 +6201,31 @@ CONTAINS 0.0_dp, m_tmp_no_1, & retain_sparsity=.TRUE.) - !!! slower way of taking the norm into account - !!CALL dbcsr_copy(m_tmp_no_1,m_tmp_no_2) - !!CALL dbcsr_multiply("N","N",1.0_dp,& - !! m_tmp_no_2,& - !! m_sig_sqrti_ii,& - !! 0.0_dp,m_tmp_no_1,& - !! retain_sparsity=.TRUE.,& - !! ) - !! - !!! get [tr(T).G]_ii - !!CALL dbcsr_copy(m_tmp_oo_1,m_sig_sqrti_ii) - !!CALL dbcsr_multiply("T","N",1.0_dp,& - !! m_t,& - !! m_tmp_no_2,& - !! 0.0_dp,m_tmp_oo_1,& - !! retain_sparsity=.TRUE.,& - !! ) - !!CALL dbcsr_get_info(m_sig_sqrti_ii, nfullrows_total=dim0 ) - !!ALLOCATE(tg_diagonal(dim0)) - !!CALL dbcsr_get_diag(m_tmp_oo_1,tg_diagonal) - !!CALL dbcsr_set(m_tmp_oo_1,0.0_dp) - !!CALL dbcsr_set_diag(m_tmp_oo_1,tg_diagonal) - !!DEALLOCATE(tg_diagonal) - !! - !!CALL dbcsr_multiply("N","N",1.0_dp,& - !! m_sig_sqrti_ii,& - !! m_tmp_oo_1,& - !! 0.0_dp,m_tmp_oo_2,& - !! filter_eps=eps_filter,& - !! ) - !!CALL dbcsr_multiply("N","N",-1.0_dp,& - !! m_ST,& - !! m_tmp_oo_2,& - !! 1.0_dp,m_tmp_no_1,& - !! retain_sparsity=.TRUE.,& - !! ) + ! get [tr(T).G]_ii + CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii) + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_t, & + m_tmp_no_2, & + 0.0_dp, m_tmp_oo_1, & + retain_sparsity=.TRUE.) + + CALL dbcsr_get_info(m_sig_sqrti_ii, nfullrows_total=dim0) + ALLOCATE (tg_diagonal(dim0)) + CALL dbcsr_get_diag(m_tmp_oo_1, tg_diagonal) + CALL dbcsr_set(m_tmp_oo_1, 0.0_dp) + CALL dbcsr_set_diag(m_tmp_oo_1, tg_diagonal) + DEALLOCATE (tg_diagonal) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_sig_sqrti_ii, & + m_tmp_oo_1, & + 0.0_dp, m_tmp_oo_2, & + filter_eps=eps_filter) + CALL dbcsr_multiply("N", "N", -1.0_dp, & + m_ST, & + m_tmp_oo_2, & + 1.0_dp, m_tmp_no_1, & + retain_sparsity=.TRUE.) ELSE @@ -5037,6 +6329,10 @@ CONTAINS CALL dbcsr_release(m_tmp_no_3) CALL dbcsr_release(m_tmp_oo_1) CALL dbcsr_release(m_tmp_oo_2) + CALL dbcsr_release(tempNOcc1) + CALL dbcsr_release(tempOccOcc1) + CALL dbcsr_release(temp1) + CALL dbcsr_release(temp2) CALL timestop(handle) @@ -5052,8 +6348,8 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE print_mathematica_matrix(matrix, filename) - TYPE(dbcsr_type), INTENT(INOUT) :: matrix - CHARACTER(len=*) :: filename + TYPE(dbcsr_type), INTENT(IN) :: matrix + CHARACTER(len=*), INTENT(IN) :: filename CHARACTER(len=*), PARAMETER :: routineN = 'print_mathematica_matrix', & routineP = moduleN//':'//routineN @@ -5151,6 +6447,327 @@ CONTAINS END SUBROUTINE print_mathematica_matrix +! ***************************************************************************** +!> \brief Compute the objective functional of NLMOs +!> \param localization_obj_function_ispin ... +!> \param penalty_func_ispin ... +!> \param penalty_vol_prefactor ... +!> \param overlap_determinant ... +!> \param m_sigma ... +!> \param nocc ... +!> \param m_B0 ... +!> \param m_theta_normalized ... +!> \param template_matrix_mo ... +!> \param weights ... +!> \param m_S0 ... +!> \param just_started ... +!> \param penalty_amplitude ... +!> \param eps_filter ... +!> \par History +!> 2020.01 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE compute_obj_nlmos(localization_obj_function_ispin, penalty_func_ispin, & + penalty_vol_prefactor, overlap_determinant, m_sigma, nocc, m_B0, & + m_theta_normalized, template_matrix_mo, weights, m_S0, just_started, & + penalty_amplitude, eps_filter) + + REAL(KIND=dp), INTENT(INOUT) :: localization_obj_function_ispin, penalty_func_ispin, & + penalty_vol_prefactor, overlap_determinant + TYPE(dbcsr_type), INTENT(INOUT) :: m_sigma + INTEGER, INTENT(IN) :: nocc + TYPE(dbcsr_type), DIMENSION(:, :), INTENT(IN) :: m_B0 + TYPE(dbcsr_type), INTENT(IN) :: m_theta_normalized, template_matrix_mo + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: weights + TYPE(dbcsr_type), INTENT(IN) :: m_S0 + LOGICAL, INTENT(IN) :: just_started + REAL(KIND=dp), INTENT(IN) :: penalty_amplitude, eps_filter + + CHARACTER(len=*), PARAMETER :: routineN = 'compute_obj_nlmos', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, idim0, ielem, para_group, reim + REAL(KIND=dp) :: det1, fval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: reim_diag, z2 + TYPE(dbcsr_type) :: tempNOcc1, tempOccOcc1, tempOccOcc2 + + CALL timeset(routineN, handle) + + CALL dbcsr_create(tempNOcc1, & + template=template_matrix_mo, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc1, & + template=m_theta_normalized, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(tempOccOcc2, & + template=m_theta_normalized, & + matrix_type=dbcsr_type_no_symmetry) + + localization_obj_function_ispin = 0.0_dp + penalty_func_ispin = 0.0_dp + ALLOCATE (z2(nocc)) + ALLOCATE (reim_diag(nocc)) + + CALL dbcsr_get_info(tempOccOcc2, group=para_group) + + DO idim0 = 1, SIZE(m_B0, 2) ! this loop is over miller ind + + z2(:) = 0.0_dp + + DO reim = 1, SIZE(m_B0, 1) ! this loop is over Re/Im + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_B0(reim, idim0), & + m_theta_normalized, & + 0.0_dp, tempOccOcc1, & + filter_eps=eps_filter) + CALL dbcsr_set(tempOccOcc2, 0.0_dp) + CALL dbcsr_add_on_diag(tempOccOcc2, 1.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_theta_normalized, & + tempOccOcc1, & + 0.0_dp, tempOccOcc2, & + retain_sparsity=.TRUE.) + + reim_diag = 0.0_dp + CALL dbcsr_get_diag(tempOccOcc2, reim_diag) + CALL mp_sum(reim_diag, para_group) + z2(:) = z2(:) + reim_diag(:)*reim_diag(:) + + ENDDO + + DO ielem = 1, nocc + SELECT CASE (2) ! allows for selection of different spread functionals + CASE (1) ! functional = -W_I * log( |z_I|^2 ) + fval = -weights(idim0)*LOG(ABS(z2(ielem))) + CASE (2) ! functional = W_I * ( 1 - |z_I|^2 ) + fval = weights(idim0) - weights(idim0)*ABS(z2(ielem)) + CASE (3) ! functional = W_I * ( 1 - |z_I| ) + fval = weights(idim0) - weights(idim0)*SQRT(ABS(z2(ielem))) + END SELECT + localization_obj_function_ispin = localization_obj_function_ispin + fval + ENDDO + + ENDDO ! end loop over idim0 + + DEALLOCATE (z2) + DEALLOCATE (reim_diag) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_S0, & + m_theta_normalized, & + 0.0_dp, tempOccOcc1, & + filter_eps=eps_filter) + ! compute current sigma + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_theta_normalized, & + tempOccOcc1, & + 0.0_dp, m_sigma, & + filter_eps=eps_filter) + + CALL determinant(m_sigma, det1, & + eps_filter) + ! save the current determinant + overlap_determinant = det1 + + IF (just_started .AND. penalty_amplitude .LT. 0.0_dp) THEN + penalty_vol_prefactor = -(-penalty_amplitude)*localization_obj_function_ispin + ENDIF + penalty_func_ispin = penalty_func_ispin + penalty_vol_prefactor*LOG(det1) + + CALL dbcsr_release(tempNOcc1) + CALL dbcsr_release(tempOccOcc1) + CALL dbcsr_release(tempOccOcc2) + + CALL timestop(handle) + + END SUBROUTINE compute_obj_nlmos + +! ***************************************************************************** +!> \brief Compute the gradient wrt the main variable +!> \param m_grad_out ... +!> \param m_B0 ... +!> \param weights ... +!> \param m_S0 ... +!> \param m_theta_normalized ... +!> \param m_siginv ... +!> \param m_sig_sqrti_ii ... +!> \param penalty_vol_prefactor ... +!> \param eps_filter ... +!> \param suggested_vol_penalty ... +!> \par History +!> 2018.10 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE compute_gradient_nlmos(m_grad_out, m_B0, weights, & + m_S0, m_theta_normalized, m_siginv, m_sig_sqrti_ii, & + penalty_vol_prefactor, eps_filter, suggested_vol_penalty) + + TYPE(dbcsr_type), INTENT(INOUT) :: m_grad_out + TYPE(dbcsr_type), DIMENSION(:, :), INTENT(IN) :: m_B0 + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: weights + TYPE(dbcsr_type), INTENT(IN) :: m_S0, m_theta_normalized, m_siginv, & + m_sig_sqrti_ii + REAL(KIND=dp), INTENT(IN) :: penalty_vol_prefactor, eps_filter + REAL(KIND=dp), INTENT(INOUT) :: suggested_vol_penalty + + CHARACTER(len=*), PARAMETER :: routineN = 'compute_gradient_nlmos', & + routineP = moduleN//':'//routineN + + INTEGER :: dim0, handle, idim0, reim + REAL(KIND=dp) :: norm_loc, norm_vol + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tg_diagonal, z2 + TYPE(dbcsr_type) :: m_temp_oo_1, m_temp_oo_2, m_temp_oo_3, & + m_temp_oo_4 + + CALL timeset(routineN, handle) + + CALL dbcsr_create(m_temp_oo_1, & + template=m_theta_normalized, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_temp_oo_2, & + template=m_theta_normalized, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_temp_oo_3, & + template=m_theta_normalized, & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_temp_oo_4, & + template=m_theta_normalized, & + matrix_type=dbcsr_type_no_symmetry) + + CALL dbcsr_get_info(m_siginv, nfullrows_total=dim0) + ALLOCATE (tg_diagonal(dim0)) + ALLOCATE (z2(dim0)) + CALL dbcsr_set(m_temp_oo_1, 0.0_dp) ! accumulate the gradient wrt a_norm here + + ! do d_Omega/d_a_normalized first + DO idim0 = 1, SIZE(m_B0, 2) ! this loop is over miller ind + + z2(:) = 0.0_dp + CALL dbcsr_set(m_temp_oo_2, 0.0_dp) ! accumulate index gradient here + DO reim = 1, SIZE(m_B0, 1) ! this loop is over Re/Im + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_B0(reim, idim0), & + m_theta_normalized, & + 0.0_dp, m_temp_oo_3, & + filter_eps=eps_filter) + + ! result contain Re/Im part of Z for the current Miller index + ! warning - save time by computing only the diagonal elements + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_theta_normalized, & + m_temp_oo_3, & + 0.0_dp, m_temp_oo_4, & + filter_eps=eps_filter) + + tg_diagonal(:) = 0.0_dp + CALL dbcsr_get_diag(m_temp_oo_4, tg_diagonal) + CALL dbcsr_set(m_temp_oo_4, 0.0_dp) + CALL dbcsr_set_diag(m_temp_oo_4, tg_diagonal) + !CALL mp_sum(tg_diagonal, para_group) + z2(:) = z2(:) + tg_diagonal(:)*tg_diagonal(:) + + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_temp_oo_3, & + m_temp_oo_4, & + 1.0_dp, m_temp_oo_2, & + filter_eps=eps_filter) + + ENDDO + + ! TODO: because some elements are zeros on some MPI tasks the + ! gradient evaluation will fail for CASE 1 and 3 + SELECT CASE (2) ! allows for selection of different spread functionals + CASE (1) ! functional = -W_I * log( |z_I|^2 ) + z2(:) = -weights(idim0)/z2(:) + CASE (2) ! functional = W_I * ( 1 - |z_I|^2 ) + z2(:) = -weights(idim0) + CASE (3) ! functional = W_I * ( 1 - |z_I| ) + z2(:) = -weights(idim0)/(2*SQRT(z2(:))) + END SELECT + CALL dbcsr_set(m_temp_oo_3, 0.0_dp) + CALL dbcsr_set_diag(m_temp_oo_3, z2) + ! TODO: print this matrix to make sure its block structure is fine + ! and there are no unecessary elements + + CALL dbcsr_multiply("N", "N", 4.0_dp, & + m_temp_oo_2, & + m_temp_oo_3, & + 1.0_dp, m_temp_oo_1, & + filter_eps=eps_filter) + + ENDDO ! end loop over idim0 + DEALLOCATE (z2) + + ! sigma0.a_norm is necessary for the volume penalty and normalization + CALL dbcsr_multiply("N", "N", & + 1.0_dp, & + m_S0, & + m_theta_normalized, & + 0.0_dp, m_temp_oo_2, & + filter_eps=eps_filter) + + ! add gradient of the penalty functional log[det(sigma)] + ! G = 2*prefactor*sigma0.a_norm.sigma_inv + CALL dbcsr_multiply("N", "N", & + 1.0_dp, & + m_temp_oo_2, & + m_siginv, & + 0.0_dp, m_temp_oo_3, & + filter_eps=eps_filter) + CALL dbcsr_norm(m_temp_oo_3, & + dbcsr_norm_maxabsnorm, norm_scalar=norm_vol) + CALL dbcsr_norm(m_temp_oo_1, & + dbcsr_norm_maxabsnorm, norm_scalar=norm_loc) + suggested_vol_penalty = norm_loc/norm_vol + CALL dbcsr_add(m_temp_oo_1, m_temp_oo_3, & + 1.0_dp, 2.0_dp*penalty_vol_prefactor) + + ! take into account the factor from the normalization constraint + ! G = ( G - sigma0.a_norm.[tr(a_norm).G]_ii ) . [sig_sqrti]_ii + ! 1. get G.[sig_sqrti]_ii + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_temp_oo_1, & + m_sig_sqrti_ii, & + 0.0_dp, m_grad_out, & + filter_eps=eps_filter) + + ! 2. get [tr(a_norm).G]_ii + ! it is possible to save time by computing only the diagonal elements + CALL dbcsr_multiply("T", "N", 1.0_dp, & + m_theta_normalized, & + m_temp_oo_1, & + 0.0_dp, m_temp_oo_3, & + filter_eps=eps_filter) + CALL dbcsr_get_diag(m_temp_oo_3, tg_diagonal) + CALL dbcsr_set(m_temp_oo_3, 0.0_dp) + CALL dbcsr_set_diag(m_temp_oo_3, tg_diagonal) + + ! 3. [X]_ii . [sig_sqrti]_ii + ! it is possible to save time by computing only the diagonal elements + CALL dbcsr_multiply("N", "N", 1.0_dp, & + m_sig_sqrti_ii, & + m_temp_oo_3, & + 0.0_dp, m_temp_oo_1, & + filter_eps=eps_filter) + ! 4. (sigma0*a_norm) .[X]_ii + CALL dbcsr_multiply("N", "N", -1.0_dp, & + m_temp_oo_2, & + m_temp_oo_1, & + 1.0_dp, m_grad_out, & + filter_eps=eps_filter) + + DEALLOCATE (tg_diagonal) + CALL dbcsr_release(m_temp_oo_1) + CALL dbcsr_release(m_temp_oo_2) + CALL dbcsr_release(m_temp_oo_3) + CALL dbcsr_release(m_temp_oo_4) + + CALL timestop(handle) + + END SUBROUTINE compute_gradient_nlmos + ! ***************************************************************************** !> \brief Compute MO coeffs from the main optimized variable (e.g. Theta, X) !> \param m_var_in ... @@ -5361,9 +6978,8 @@ CONTAINS TYPE(domain_submatrix_type), DIMENSION(:), & INTENT(INOUT) :: domain_prec_out TYPE(dbcsr_type), INTENT(INOUT) :: m_prec_out, m_ks, m_s - TYPE(dbcsr_type), INTENT(IN) :: m_siginv - TYPE(dbcsr_type), INTENT(INOUT) :: m_quench_t - TYPE(dbcsr_type), INTENT(IN) :: m_FTsiginv, m_siginvTFTsiginv, m_ST + TYPE(dbcsr_type), INTENT(IN) :: m_siginv, m_quench_t, m_FTsiginv, & + m_siginvTFTsiginv, m_ST TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: m_STsiginv_out, m_s_vv_out, m_f_vv_out TYPE(cp_para_env_type), POINTER :: para_env TYPE(cp_blacs_env_type), POINTER :: blacs_env @@ -5835,7 +7451,7 @@ CONTAINS penalty_occ_vol, normalize_orbitals, penalty_occ_vol_prefactor, & penalty_occ_vol_pf2, special_case) - TYPE(optimizer_options_type) :: optimizer + TYPE(optimizer_options_type), INTENT(IN) :: optimizer TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_grad TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_delta, m_s, m_ks, m_siginv, m_quench_t TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_FTsiginv, m_siginvTFTsiginv, m_ST, & @@ -6390,10 +8006,8 @@ CONTAINS normalize_orbitals, penalty_occ_vol_prefactor, eps_filter, path_num) TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_x_in, m_x_out, m_ks, m_s - TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_siginv - TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_quench_t - TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_FTsiginv, m_siginvTFTsiginv, m_ST, & - m_STsiginv + TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_siginv, m_quench_t, m_FTsiginv, & + m_siginvTFTsiginv, m_ST, m_STsiginv TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_s_vv, m_ks_vv, m_g_full TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: m_t, m_sig_sqrti_ii LOGICAL, INTENT(IN) :: penalty_occ_vol, normalize_orbitals @@ -7449,8 +9063,8 @@ CONTAINS special_case) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(almo_scf_env_type) :: almo_scf_env - TYPE(optimizer_options_type) :: optimizer + TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env + TYPE(optimizer_options_type), INTENT(IN) :: optimizer TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: quench_t, matrix_t_in, matrix_t_out LOGICAL, INTENT(IN) :: perturbation_only INTEGER, INTENT(IN), OPTIONAL :: special_case @@ -7537,10 +9151,11 @@ CONTAINS ENDIF ! penalty amplitude adjusts the strenght of volume conservation - penalty_occ_vol = (almo_scf_env%penalty%occ_vol_method .NE. almo_occ_vol_penalty_none .AND. & - my_special_case .EQ. xalmo_case_fully_deloc) + penalty_occ_vol = .FALSE. + !(almo_scf_env%penalty%occ_vol_method .NE. almo_occ_vol_penalty_none .AND. & + ! my_special_case .EQ. xalmo_case_fully_deloc) normalize_orbitals = penalty_occ_vol - penalty_amplitude = almo_scf_env%penalty%occ_vol_coeff + penalty_amplitude = 0.0_dp !almo_scf_env%penalty%occ_vol_coeff ALLOCATE (penalty_occ_vol_g_prefactor(nspins)) ALLOCATE (penalty_occ_vol_h_prefactor(nspins)) penalty_occ_vol_g_prefactor(:) = 0.0_dp @@ -8692,7 +10307,7 @@ CONTAINS spin_factor = 1.0_dp ENDIF - penalty_amplitude = almo_scf_env%penalty%occ_vol_coeff + penalty_amplitude = 0.0_dp !almo_scf_env%penalty%occ_vol_coeff ALLOCATE (nocc(nspins)) DO ispin = 1, nspins diff --git a/src/almo_scf_types.F b/src/almo_scf_types.F index e615f0571f..bf93ba634b 100644 --- a/src/almo_scf_types.F +++ b/src/almo_scf_types.F @@ -46,8 +46,11 @@ MODULE almo_scf_types ! methods that add penalty terms to the energy functional TYPE penalty_type - REAL(KIND=dp) :: occ_vol_coeff - INTEGER :: occ_vol_method + REAL(KIND=dp) :: final_determinant, penalty_strength, & + determinant_tolerance, penalty_strength_dec_factor, & + compactification_filter_start + INTEGER :: operator_type + LOGICAL :: virtual_nlmos END TYPE penalty_type @@ -74,6 +77,7 @@ MODULE almo_scf_types neglect_threshold INTEGER :: optimizer_type ! diis, pcg, etc. + TYPE(penalty_type) :: opt_penalty INTEGER :: preconditioner, & ! preconditioner type conjugator, & ! conjugator type @@ -207,7 +211,7 @@ MODULE almo_scf_types LOGICAL :: s_sqrt_done REAL(KIND=dp) :: almo_scf_energy LOGICAL :: orthogonal_basis, fixed_mu - LOGICAL :: return_orthogonalized_mos + LOGICAL :: return_orthogonalized_mos, construct_nlmos !! Smearing control !! smear flag allow to retrieve eigenvalues in almo_scf with diag algorithm and create occupation-scaled ALMO orbitals @@ -239,6 +243,9 @@ MODULE almo_scf_types ! mo overlap inversion algorithm INTEGER :: sigma_inv_algorithm + ! Determinant of the ALMO overlap matrix + REAL(KIND=dp) :: overlap_determinant + ! ALMO SCF delocalization control LOGICAL :: perturbative_delocalization INTEGER :: quencher_radius_type @@ -354,7 +361,6 @@ MODULE almo_scf_types ! Options for various subsection options collected neatly TYPE(almo_analysis_type) :: almo_analysis - TYPE(penalty_type) :: penalty ! Options for various optimizers collected neatly TYPE(optimizer_options_type) :: opt_block_diag_diis @@ -362,6 +368,7 @@ MODULE almo_scf_types TYPE(optimizer_options_type) :: opt_xalmo_diis TYPE(optimizer_options_type) :: opt_xalmo_pcg TYPE(optimizer_options_type) :: opt_xalmo_trustr + TYPE(optimizer_options_type) :: opt_nlmo_pcg TYPE(optimizer_options_type) :: opt_xalmo_newton_pcg_solver TYPE(optimizer_options_type) :: opt_k_pcg diff --git a/src/input_constants.F b/src/input_constants.F index 743f04b9e2..51d3f2cd8e 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -917,8 +917,9 @@ MODULE input_constants almo_scf_trustr = 4, & almo_scf_skip = 0 - INTEGER, PARAMETER, PUBLIC :: almo_occ_vol_penalty_none = 0, & - almo_occ_vol_penalty_lndet = 1 + INTEGER, PARAMETER, PUBLIC :: penalty_type_none = 0, & + penalty_type_lndet = 1, & + penalty_type_nlmo = 2 ! optimizer parameters INTEGER, PARAMETER, PUBLIC :: cg_zero = 0, & diff --git a/src/input_cp2k_almo.F b/src/input_cp2k_almo.F index 2f3542de60..93d0fcd7c2 100644 --- a/src/input_cp2k_almo.F +++ b/src/input_cp2k_almo.F @@ -17,13 +17,13 @@ MODULE input_cp2k_almo USE input_constants, ONLY: & almo_deloc_none, almo_deloc_scf, almo_deloc_x, almo_deloc_x_then_scf, & almo_deloc_xalmo_1diag, almo_deloc_xalmo_scf, almo_deloc_xalmo_x, almo_frz_crystal, & - almo_frz_none, almo_occ_vol_penalty_lndet, almo_occ_vol_penalty_none, almo_scf_diag, & - almo_scf_pcg, almo_scf_skip, almo_scf_trustr, atomic_guess, cg_dai_yuan, cg_fletcher, & - cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, cg_liu_storey, cg_polak_ribiere, & - cg_zero, molecular_guess, optimizer_diis, optimizer_lin_eq_pcg, optimizer_pcg, & - optimizer_trustr, spd_inversion_dense_cholesky, spd_inversion_ls_hotelling, & - spd_inversion_ls_taylor, trustr_cauchy, trustr_dogleg, trustr_steihaug, xalmo_prec_domain, & - xalmo_prec_full, xalmo_prec_zero, xalmo_trial_r0_out, xalmo_trial_simplex + almo_frz_none, almo_scf_diag, almo_scf_pcg, almo_scf_skip, almo_scf_trustr, atomic_guess, & + cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, & + cg_liu_storey, cg_polak_ribiere, cg_zero, molecular_guess, op_loc_berry, op_loc_pipek, & + optimizer_diis, optimizer_lin_eq_pcg, optimizer_pcg, optimizer_trustr, & + spd_inversion_dense_cholesky, spd_inversion_ls_hotelling, spd_inversion_ls_taylor, & + trustr_cauchy, trustr_dogleg, trustr_steihaug, xalmo_prec_domain, xalmo_prec_full, & + xalmo_prec_zero, xalmo_trial_r0_out, xalmo_trial_simplex USE input_keyword_types, ONLY: keyword_create,& keyword_release,& keyword_type @@ -46,6 +46,7 @@ MODULE input_cp2k_almo INTEGER, PARAMETER, PRIVATE :: optimizer_xalmo_pcg = 3 INTEGER, PARAMETER, PRIVATE :: optimizer_xalmo_trustr = 4 INTEGER, PARAMETER, PRIVATE :: optimizer_newton_pcg_solver = 5 + INTEGER, PARAMETER, PRIVATE :: optimizer_nlmo_pcg = 6 PUBLIC :: create_almo_scf_section @@ -210,6 +211,12 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="CONSTRUCT_NLMOS", & + description="Turns on post-SCF construction of NLMOs", & + usage="CONSTRUCT_NLMOS .TRUE.", default_l_val=.FALSE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + !CALL keyword_create(keyword, __LOCATION__, name="DOMAIN_LAYOUT_MOS",& ! description="Each electron in the system is constrained to its own delocalization domain."//& ! " This keyword creates groups of electrons that share the same domain.",& @@ -533,7 +540,7 @@ CONTAINS CALL section_release(subsection) NULLIFY (subsection) - CALL create_penalty_section(subsection) + CALL create_optimizer_section(subsection, optimizer_nlmo_pcg) CALL section_add_subsection(section, subsection) CALL section_release(subsection) @@ -584,12 +591,21 @@ CONTAINS CASE (optimizer_block_diagonal_pcg) CALL section_create(section, __LOCATION__, name="ALMO_OPTIMIZER_PCG", & description="Controls the PCG optimization of block-diagonal ALMOs.", & - n_keywords=9, n_subsections=0, repeats=.FALSE.) + n_keywords=9, n_subsections=1, repeats=.FALSE.) optimizer_type = optimizer_pcg + CASE (optimizer_nlmo_pcg) + CALL section_create(section, __LOCATION__, name="NLMO_OPTIMIZER_PCG", & + description="Controls the PCG optimization of nonorthogonal localized MOs.", & + n_keywords=9, n_subsections=1, repeats=.FALSE.) + optimizer_type = optimizer_pcg + NULLIFY (subsection) + CALL create_penalty_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) CASE (optimizer_xalmo_pcg) CALL section_create(section, __LOCATION__, name="XALMO_OPTIMIZER_PCG", & description="Controls the PCG optimization of extended ALMOs.", & - n_keywords=10, n_subsections=1, repeats=.FALSE.) + n_keywords=10, n_subsections=2, repeats=.FALSE.) NULLIFY (subsection) CALL create_optimizer_section(subsection, optimizer_newton_pcg_solver) CALL section_add_subsection(section, subsection) @@ -848,25 +864,54 @@ CONTAINS CALL section_create(section, __LOCATION__, name="PENALTY", & description="Add penalty terms to the energy functional.", & - n_keywords=2, n_subsections=0, repeats=.FALSE.) + n_keywords=3, n_subsections=0, repeats=.FALSE.) NULLIFY (keyword) CALL keyword_create( & - keyword, __LOCATION__, name="OCCUPIED_VOLUME_PENALTY_METHOD", & - description="Penalty that prevents nonorthogonal orbitals from becoming linear dependent.", & - usage="OCCUPIED_VOLUME_PENALTY_METHOD LNDET", & - default_i_val=almo_occ_vol_penalty_none, & - enum_c_vals=s2a("NONE", "LNDET"), & - enum_desc=s2a("Do not use penalties", & - "Use -coeff*ln(det(MO-overlap)) term."), & - enum_i_vals=(/almo_occ_vol_penalty_none, almo_occ_vol_penalty_lndet/)) + keyword, __LOCATION__, name="OPERATOR", & + description="Type of opertator which defines the spread functional", & + usage="OPERATOR PIPEK", & + enum_c_vals=s2a("BERRY", "PIPEK"), & + enum_i_vals=(/op_loc_berry, op_loc_pipek/), & + default_i_val=op_loc_berry) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="OCCUPIED_VOLUME_PENALTY_COEFF", & - description="Multiplication factor that determines the strength of the penalty term.", & - usage="OCCUPIED_VOLUME_PENALTY_COEFF 10.0", default_r_val=1.0_dp) + CALL keyword_create(keyword, __LOCATION__, name="PENALTY_STRENGTH", & + description="Strength of the orthogonalization penalty", & + usage="PENALTY_STRENGTH 1.1", default_r_val=1.1_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="PENALTY_STRENGTH_DECREASE_FACTOR", & + description="Factor that decreases the strength of the orthogonalization penalty.", & + usage="PENALTY_STRENGTH_DECREASE_FACTOR 1.1", default_r_val=1.1_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="DETERMINANT_TOLERANCE", & + description="Stop the optimization of the penalty strength if the determinant of the overlap "// & + "changes less than this tolerance threshold.", & + usage="DETERMINANT_TOLERANCE 1.0E-4", default_r_val=1.0E-3_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="FINAL_DETERMINANT", & + description="The final determinant that obtained after optimization.", & + usage="FINAL_DETERMINANT 0.1", default_r_val=0.1_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="COMPACTIFICATION_FILTER_START", & + description="Set orbital coefficients with absolute value smaller than this value to zero.", & + usage="COMPACTIFICATION_FILTER_START 1.e-6", default_r_val=-1.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="VIRTUAL_NLMOS", & + description="Localize virtual oribtals", & + usage="VIRTUAL_NLMOS .TRUE.", default_l_val=.FALSE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/iterate_matrix.F b/src/iterate_matrix.F index 0f44b91e91..251f5439d6 100644 --- a/src/iterate_matrix.F +++ b/src/iterate_matrix.F @@ -28,6 +28,7 @@ MODULE iterate_matrix m_walltime USE mathconstants, ONLY: ifac USE mathlib, ONLY: abnormal_value + USE message_passing, ONLY: mp_sum #include "./base/base_uses.f90" IMPLICIT NONE @@ -70,8 +71,9 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'determinant', routineP = moduleN//':'//routineN - INTEGER :: handle, i, max_iter_lanczos, nsize, & - order_lanczos, sign_iter, unit_nr + INTEGER :: group, handle, i, max_iter_lanczos, & + nsize, order_lanczos, sign_iter, & + unit_nr INTEGER(KIND=int_8) :: flop1 INTEGER, SAVE :: recursion_depth = 0 REAL(KIND=dp) :: det0, eps_lanczos, frobnorm, maxnorm, & @@ -102,9 +104,10 @@ CONTAINS matrix_type=dbcsr_type_no_symmetry) ! compute the product of the diagonal elements - CALL dbcsr_get_info(matrix, nfullrows_total=nsize) + CALL dbcsr_get_info(matrix, nfullrows_total=nsize, group=group) ALLOCATE (diagonal(nsize)) CALL dbcsr_get_diag(matrix, diagonal) + CALL mp_sum(diagonal, group) det = PRODUCT(diagonal) ! create diagonal SQRTI matrix diff --git a/src/qs_loc_methods.F b/src/qs_loc_methods.F index 3e24d71cb3..eb5f2707a5 100644 --- a/src/qs_loc_methods.F +++ b/src/qs_loc_methods.F @@ -302,6 +302,7 @@ CONTAINS NULLIFY (particle_set) ! get rows and cols of the input CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc) + ! replicate the input kind of matrix CALL cp_fm_create(opvec, vectors%matrix_struct) CALL cp_fm_set_all(opvec, 0.0_dp) @@ -329,18 +330,22 @@ CONTAINS CALL cp_fm_get_info(zij_fm_set(iatom, 1)%matrix, ncol_global=ldz) isgf = first_sgf(iatom) ncol = nsgf(iatom) + ! multiply fmxfm, using only part of the ao : Ct x S CALL cp_gemm('N', 'N', nao, nmoloc, nao, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, & a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1) + CALL cp_gemm('T', 'N', nmoloc, nmoloc, ncol, 0.5_dp, vectors, opvec, & 0.0_dp, zij_fm_set(iatom, 1)%matrix, & a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf) CALL cp_gemm('N', 'N', nao, nmoloc, ncol, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, & - a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf) + a_first_col=isgf, a_first_row=1, b_first_col=1, b_first_row=isgf) + CALL cp_gemm('T', 'N', nmoloc, nmoloc, nao, 0.5_dp, vectors, opvec, & 1.0_dp, zij_fm_set(iatom, 1)%matrix, & a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1) + END DO ! iatom ! And now perform the optimization and rotate the orbitals @@ -646,9 +651,12 @@ CONTAINS END IF DO ispin = s_spin, l_spin + CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) loc_method = localized_wfn_control%localization_method + SELECT CASE (localized_wfn_control%operator_type) + CASE (op_loc_berry) ! Here we allocate op_fm_set with the RIGHT size for uks NULLIFY (tmp_fm_struct, mo_coeff) @@ -691,6 +699,7 @@ CONTAINS CPABORT("Boys localization not implemented") CASE (op_loc_pipek) + CALL optimize_loc_pipek(qs_env, loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, & op_fm_set, ispin, print_loc_section) diff --git a/src/qs_loc_utils.F b/src/qs_loc_utils.F index ea6b3d71c7..db74401685 100644 --- a/src/qs_loc_utils.F +++ b/src/qs_loc_utils.F @@ -99,7 +99,7 @@ MODULE qs_loc_utils ! *** Public *** PUBLIC :: jacobi_rotation_pipek, qs_loc_env_init, loc_write_restart, & - retain_history, qs_loc_init, & + retain_history, qs_loc_init, compute_berry_operator, & set_loc_centers, set_loc_wfn_lists, qs_loc_control_init CONTAINS @@ -278,7 +278,9 @@ CONTAINS ct = COS(theta) st = SIN(theta) + CALL rotate_rmat_real(istate, jstate, st, ct, rmat) + CALL rotate_zij_real(istate, jstate, st, ct, zij_fm_set) END DO END DO @@ -649,21 +651,13 @@ CONTAINS END SUBROUTINE qs_loc_env_init ! ************************************************************************************************** -!> \brief Computes the Berry operator for periodic systems -!> used to define the spread of the MOS -!> Here the matrix elements of the type and -!> are computed, where mu and nu are the contracted basis functions. -!> Namely the Berry operator is exp(ikr) -!> k is defined somewhere -!> the pair lists are exploited and sparse matrixes are constructed +!> \brief A wrapper to compute the Berry operator for periodic systems !> \param qs_loc_env new environment for the localization calculations !> \param qs_env the qs_env in which the qs_env lives !> \par History !> 04.2005 created [MI] +!> 04.2018 modified [RZK, ZL] !> \author MI -!> \note -!> The intgrals are computed analytically using the primitives GTO -!> The contraction is performed block-wise ! ************************************************************************************************** SUBROUTINE get_berry_operator(qs_loc_env, qs_env) TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env @@ -672,9 +666,51 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'get_berry_operator', & routineP = moduleN//':'//routineN - INTEGER :: dim_op, handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, & - last_jatom, ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, & - sgfb + INTEGER :: dim_op, handle + TYPE(cell_type), POINTER :: cell + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set + + CALL timeset(routineN, handle) + + NULLIFY (cell, op_sm_set) + CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, & + cell=cell, dim_op=dim_op) + CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op) + + CALL timestop(handle) + END SUBROUTINE get_berry_operator + +! ************************************************************************************************** +!> \brief Computes the Berry operator for periodic systems +!> used to define the spread of the MOS +!> Here the matrix elements of the type and +!> are computed, where mu and nu are the contracted basis functions. +!> Namely the Berry operator is exp(ikr) +!> k is defined somewhere +!> the pair lists are exploited and sparse matrixes are constructed +!> \param qs_env the qs_env in which the qs_env lives +!> \param cell ... +!> \param op_sm_set ... +!> \param dim_op ... +!> \par History +!> 04.2005 created [MI] +!> 04.2018 wrapped old code [RZK, ZL] +!> \author MI +!> \note +!> The intgrals are computed analytically using the primitives GTO +!> The contraction is performed block-wise +! ************************************************************************************************** + SUBROUTINE compute_berry_operator(qs_env, cell, op_sm_set, dim_op) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(cell_type), POINTER :: cell + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set + INTEGER :: dim_op + + CHARACTER(len=*), PARAMETER :: routineN = 'compute_berry_operator', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, & + ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, sgfb INTEGER, DIMENSION(3) :: perd0 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & npgfb, nsgfa, nsgfb @@ -686,8 +722,6 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, & sphi_b, work, zeta, zetb TYPE(block_p_type), DIMENSION(:), POINTER :: op_cos, op_sin - TYPE(cell_type), POINTER :: cell - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b TYPE(neighbor_list_iterator_p_type), & @@ -699,9 +733,8 @@ CONTAINS TYPE(qs_kind_type), POINTER :: qs_kind CALL timeset(routineN, handle) - NULLIFY (qs_kind, qs_kind_set) - NULLIFY (cell, op_sm_set, particle_set) + NULLIFY (particle_set) NULLIFY (sab_orb) NULLIFY (cosab, sinab, work) NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb) @@ -722,9 +755,6 @@ CONTAINS ALLOCATE (work(ldwork, ldwork)) work = 0.0_dp - CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, & - cell=cell, dim_op=dim_op) - ALLOCATE (op_cos(dim_op)) ALLOCATE (op_sin(dim_op)) DO i = 1, dim_op @@ -838,7 +868,6 @@ CONTAINS ! *** Calculate the primitive overlap integrals *** DO i = 1, dim_op - kvec(1:3) = vector_k(1:3, i) cosab = 0.0_dp sinab = 0.0_dp @@ -846,7 +875,6 @@ CONTAINS la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), & rpgfb(:, jset), lb_min(jset), & ra, rb, kvec, cosab, sinab) - CALL contract_cossin(op_cos(i)%block, op_sin(i)%block, & iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, & jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, & @@ -876,7 +904,7 @@ CONTAINS DEALLOCATE (cosab, sinab, work, basis_set_list) CALL timestop(handle) - END SUBROUTINE get_berry_operator + END SUBROUTINE compute_berry_operator ! ************************************************************************************************** !> \brief ... diff --git a/tests/QS/regtest-nlmo/Si-nlmos.inp b/tests/QS/regtest-nlmo/Si-nlmos.inp new file mode 100644 index 0000000000..6185d27ad6 --- /dev/null +++ b/tests/QS/regtest-nlmo/Si-nlmos.inp @@ -0,0 +1,93 @@ +&GLOBAL + PROJECT NLMOS + RUN_TYPE ENERGY + PRINT_LEVEL LOW +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + + &MGRID + CUTOFF 200 + NGRIDS 4 + &END MGRID + + &QS + ALMO_SCF T + EPS_DEFAULT 1.0E-08 + &END QS + + &ALMO_SCF + + EPS_FILTER 1.0E-09 + ALMO_ALGORITHM SKIP + MO_OVERLAP_INV_ALG DENSE_CHOLESKY + DELOCALIZE_METHOD FULL_SCF + ALMO_SCF_GUESS ATOMIC + XALMO_TRIAL_WF SIMPLE + CONSTRUCT_NLMOS TRUE + + &XALMO_OPTIMIZER_PCG + MAX_ITER 50 + EPS_ERROR 1.0E-3 + CONJUGATOR HESTENES_STIEFEL + PRECONDITIONER DEFAULT + LIN_SEARCH_EPS_ERROR 0.1 + LIN_SEARCH_STEP_SIZE_GUESS 0.2 + MAX_ITER_OUTER_LOOP 10 + &END XALMO_OPTIMIZER_PCG + + &NLMO_OPTIMIZER_PCG + MAX_ITER 200 + EPS_ERROR 0.1 + CONJUGATOR ZERO + PRECONDITIONER DEFAULT + LIN_SEARCH_EPS_ERROR 0.5 + LIN_SEARCH_STEP_SIZE_GUESS 0.001 + MAX_ITER_OUTER_LOOP 0 + &PENALTY + OPERATOR PIPEK + PENALTY_STRENGTH 0.01 + DETERMINANT_TOLERANCE 1.0E-8 + PENALTY_STRENGTH_DECREASE_FACTOR 2.0 + FINAL_DETERMINANT 0.6 + COMPACTIFICATION_FILTER_START 1.0E-2 + VIRTUAL_NLMOS FALSE + &END PENALTY + &END NLMO_OPTIMIZER_PCG + + &END ALMO_SCF + + &XC + &XC_FUNCTIONAL BLYP + &END XC_FUNCTIONAL + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 5.430710 5.430710 5.430710 + &END CELL + &COORD + SCALED T + Si 0.000000 0.000000 0.000000 + Si 0.000000 0.500000 0.500000 + Si 0.500000 0.000000 0.500000 + Si 0.500000 0.500000 0.000000 + Si 0.250000 0.250000 0.250000 + Si 0.250000 0.750000 0.750000 + Si 0.750000 0.250000 0.750000 + Si 0.750000 0.750000 0.250000 + &END COORD + &KIND Si + BASIS_SET SZV-GTH + POTENTIAL GTH-PBE-q4 + &END KIND + &END SUBSYS + +&END FORCE_EVAL + diff --git a/tests/QS/regtest-nlmo/TEST_FILES b/tests/QS/regtest-nlmo/TEST_FILES new file mode 100644 index 0000000000..a3532474c4 --- /dev/null +++ b/tests/QS/regtest-nlmo/TEST_FILES @@ -0,0 +1,3 @@ +pipek_C6H6.inp 90 1e-04 177.2020586640 +Si-nlmos.inp 90 1e-04 118.5296581594 +#EOF diff --git a/tests/QS/regtest-nlmo/pipek_C6H6.inp b/tests/QS/regtest-nlmo/pipek_C6H6.inp new file mode 100644 index 0000000000..33dd70c83d --- /dev/null +++ b/tests/QS/regtest-nlmo/pipek_C6H6.inp @@ -0,0 +1,102 @@ +&GLOBAL + PROJECT C6H6 + RUN_TYPE ENERGY + PRINT_LEVEL LOW +&END GLOBAL +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 200 + NGRIDS 4 + &END MGRID + &QS + ALMO_SCF T + EPS_DEFAULT 1.0E-8 + &END QS + + &ALMO_SCF + + EPS_FILTER 1.0E-09 + ALMO_ALGORITHM SKIP + MO_OVERLAP_INV_ALG DENSE_CHOLESKY + DELOCALIZE_METHOD FULL_SCF + ALMO_SCF_GUESS ATOMIC + XALMO_TRIAL_WF SIMPLE + CONSTRUCT_NLMOS TRUE + RETURN_ORTHOGONALIZED_MOS FALSE + + &XALMO_OPTIMIZER_PCG + MAX_ITER 50 + EPS_ERROR 1.0E-3 + CONJUGATOR HESTENES_STIEFEL + PRECONDITIONER DEFAULT + LIN_SEARCH_EPS_ERROR 0.01 + LIN_SEARCH_STEP_SIZE_GUESS 0.2 + MAX_ITER_OUTER_LOOP 2 + &END XALMO_OPTIMIZER_PCG + + &NLMO_OPTIMIZER_PCG + MAX_ITER 20000 + EPS_ERROR 0.2 + CONJUGATOR ZERO + PRECONDITIONER DEFAULT + LIN_SEARCH_EPS_ERROR 0.1 + LIN_SEARCH_STEP_SIZE_GUESS 0.0001 + MAX_ITER_OUTER_LOOP 5 + &PENALTY + OPERATOR PIPEK + PENALTY_STRENGTH 0.01 + DETERMINANT_TOLERANCE 1.0E-8 + PENALTY_STRENGTH_DECREASE_FACTOR 2.0 + FINAL_DETERMINANT 0.5 + !COMPACTIFICATION_FILTER_START 1.0E-3 + VIRTUAL_NLMOS FALSE + &END PENALTY + &END NLMO_OPTIMIZER_PCG + + &END ALMO_SCF + + &XC + &XC_FUNCTIONAL BLYP + &END XC_FUNCTIONAL + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 9.0 9.0 5.0 + &END CELL + &TOPOLOGY + &CENTER_COORDINATES + &END + &GENERATE + CREATE_MOLECULES TRUE + &END + &END + &COORD + C 3.5560000000 4.5610000000 0.0000000000 + C 4.5060000000 3.5350000000 0.0000000000 + C 5.8680000000 3.8450000000 0.0000000000 + C 6.2820000000 5.1790000000 0.0000000000 + C 5.3320000000 6.2050000000 0.0000000000 + C 3.9700000000 5.8950000000 0.0000000000 + H 2.5000000000 4.3210000000 0.0000000000 + H 4.1850000000 2.5000000000 0.0000000000 + H 6.6040000000 3.0490000000 0.0000000000 + H 7.3390000000 5.4190000000 0.0000000000 + H 5.6530000000 7.2400000000 0.0000000000 + H 3.2340000000 6.6910000000 0.0000000000 + &END COORD + &KIND C + BASIS_SET SZV-GTH + POTENTIAL GTH-BLYP-q4 + &END KIND + &KIND H + BASIS_SET SZV-GTH + POTENTIAL GTH-BLYP-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 012ea57efe..90c451068b 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -263,3 +263,4 @@ QS/regtest-sto QS/regtest-ri-laplace-mp2-cubic libint mpiranks==2||mpiranks==4||mpiranks==6||mpiranks==16||mpiranks==24 QS/regtest-dm-ls-scf-4 QS/regtest-almo-trustr +QS/regtest-nlmo diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 9cd2212dd4..90e5d94888 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -1,4 +1,4 @@ -89 +90 Total energy:!3 POTENTIAL ENERGY!4 Total energy \[eV\]:!4 @@ -88,6 +88,7 @@ DIPOLE : CheckSum =!5 POLAR : CheckSum =!5 XAS excitation energy (eV): !7 Electronic density on regular grids:!7 +Final localization: !3 # # these are the tests the can be selected for regtesting. # do regtest will grep for test_grep (first column) and look if the numeric value