diff --git a/src/almo_scf.F b/src/almo_scf.F index 539afc180ae..6e4c50ca1d2 100644 --- a/src/almo_scf.F +++ b/src/almo_scf.F @@ -15,8 +15,9 @@ MODULE almo_scf almo_scf_t_rescaling,& almo_scf_t_to_proj,& distribute_domains,& - fill_matrix_with_ones,& - orthogonalize_mos + orthogonalize_mos,& + copy_virt2occ,& + fill_matrix_with_ones USE almo_scf_optimizer, ONLY: almo_scf_block_diagonal,& almo_scf_xalmo_eigensolver,& almo_scf_xalmo_pcg,& @@ -57,7 +58,8 @@ MODULE almo_scf USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: & dbcsr_add, dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, & - dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_info, dbcsr_init_random, & + dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, & + dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_init_random, dbcsr_print, & dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, & dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_nblkrows_total, & dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, & @@ -208,6 +210,7 @@ SUBROUTINE almo_scf_init(qs_env, almo_scf_env, calc_forces) 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_nlmo_trustr%optimizer_type = optimizer_trustr almo_scf_env%opt_xalmo_newton_pcg_solver%optimizer_type = optimizer_lin_eq_pcg ! get info from the qs_env @@ -524,6 +527,7 @@ SUBROUTINE almo_scf_initial_guess(qs_env, almo_scf_env) routineP = moduleN//':'//routineN CHARACTER(LEN=default_path_length) :: file_name, project_name + TYPE(dbcsr_distribution_type) :: dist INTEGER :: handle, iaspc, ispin, istore, naspc, & nspins, unit_nr INTEGER, DIMENSION(2) :: nelectron_spin @@ -532,7 +536,6 @@ SUBROUTINE almo_scf_initial_guess(qs_env, almo_scf_env) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cp_logger_type), POINTER :: logger TYPE(cp_para_env_type), POINTER :: para_env - TYPE(dbcsr_distribution_type) :: dist TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao TYPE(dft_control_type), POINTER :: dft_control TYPE(molecular_scf_guess_env_type), POINTER :: mscfg_env @@ -1680,8 +1683,10 @@ SUBROUTINE construct_nlmos(qs_env, almo_scf_env) max_iter_lanczos=almo_scf_env%max_iter_lanczos) ENDDO - CALL nlmo_optimization_entry(qs_env, almo_scf_env, & - virtuals=.FALSE.) + IF (almo_scf_env%occupied_nlmos) THEN + CALL nlmo_optimization_entry(qs_env, almo_scf_env, & + virtuals=.FALSE.) + ENDIF IF (almo_scf_env%virtual_nlmos) THEN CALL construct_virtuals(almo_scf_env) @@ -1739,31 +1744,77 @@ SUBROUTINE construct_virtuals(almo_scf_env) 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) + + CALL almo_scf_t_to_proj(t=almo_scf_env%matrix_t(ispin), & + p=almo_scf_env%matrix_p(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., & + smear=almo_scf_env%smear, & + algorithm=almo_scf_env%sigma_inv_algorithm, & + inverse_accelerator=almo_scf_env%order_lanczos, & + inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, & + 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) + 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, & + CALL dbcsr_multiply("N", "N", -1.0_dp, & + almo_scf_env%matrix_p(ispin), & tempNV1, & - almo_scf_env%matrix_t(ispin), & - 0.0_dp, tempVOcc1, & + 1.0_dp, & + almo_scf_env%matrix_v(ispin), & 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) + !! Test V.T S T = 0 + !CALL dbcsr_multiply("N", "N", 1.0_dp, & + ! almo_scf_env%matrix_s(1), & + ! almo_scf_env%matrix_t(ispin), & + ! 0.0_dp, tempNOcc, & + ! 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_multiply("T", "N", 1.0_dp, & + ! almo_scf_env%matrix_v(ispin), & + ! tempNOcc, & + ! 0.0_dp, tempVOcc1, & + ! filter_eps=almo_scf_env%eps_filter) - CALL dbcsr_add(almo_scf_env%matrix_v(ispin), tempNV1, 1.0_dp, -1.0_dp) + !CALL dbcsr_print(tempVOcc1) ! compute VxV overlap CALL dbcsr_multiply("N", "N", 1.0_dp, & @@ -1818,6 +1869,51 @@ SUBROUTINE construct_virtuals(almo_scf_env) CALL dbcsr_copy(almo_scf_env%matrix_v(ispin), tempNV1) + !! Test V.T F V = 0 + !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_print(tempVV1) + + !! Test V.T F T = 0 + !CALL dbcsr_multiply("N", "N", 1.0_dp, & + ! almo_scf_env%matrix_ks(ispin), & + ! almo_scf_env%matrix_t(ispin), & + ! 0.0_dp, tempNOcc, & + ! filter_eps=almo_scf_env%eps_filter) + + !CALL dbcsr_multiply("T", "N", 1.0_dp, & + ! almo_scf_env%matrix_v(ispin), & + ! tempNOcc, & + ! 0.0_dp, tempVOcc1, & + ! filter_eps=almo_scf_env%eps_filter) + + !CALL dbcsr_print(tempVOcc1) + + !! Test V.T S T = 0 + !CALL dbcsr_multiply("N", "N", 1.0_dp, & + ! almo_scf_env%matrix_s(1), & + ! almo_scf_env%matrix_t(ispin), & + ! 0.0_dp, tempNOcc, & + ! filter_eps=almo_scf_env%eps_filter) + + !CALL dbcsr_multiply("T", "N", 1.0_dp, & + ! almo_scf_env%matrix_v(ispin), & + ! tempNOcc, & + ! 0.0_dp, tempVOcc1, & + ! filter_eps=almo_scf_env%eps_filter) + + !CALL dbcsr_print(tempVOcc1) + CALL dbcsr_release(tempNV1) CALL dbcsr_release(tempVOcc1) CALL dbcsr_release(tempVOcc2) @@ -2051,6 +2147,15 @@ SUBROUTINE almo_scf_post(qs_env, almo_scf_env) DO ispin = 1, almo_scf_env%nspins + IF (debug_mode) THEN + !CALL dbcsr_print(almo_scf_env%matrix_t(ispin)) + !CALL dbcsr_print(almo_scf_env%matrix_v(ispin)) + CALL copy_virt2occ(m_occ=almo_scf_env%matrix_t(ispin), & + m_virt=almo_scf_env%matrix_v(ispin), & + first_virtual_of_domain = 1) + !CALL dbcsr_print(almo_scf_env%matrix_t(ispin)) + ENDIF + CALL dbcsr_create(matrix_t_processed(ispin), & template=almo_scf_env%matrix_t(ispin), & matrix_type=dbcsr_type_no_symmetry) diff --git a/src/almo_scf_env_methods.F b/src/almo_scf_env_methods.F index d156e7adb16..2ce13b5eb1c 100644 --- a/src/almo_scf_env_methods.F +++ b/src/almo_scf_env_methods.F @@ -107,8 +107,8 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env) INTEGER :: handle TYPE(section_vals_type), POINTER :: almo_analysis_section, almo_opt_diis_section, & almo_opt_pcg_section, almo_scf_section, matrix_iterate_section, nlmo_opt_pcg_section, & - nlmo_opt_trustr_section, penalty_section, xalmo_opt_newton_pcg_section, & - xalmo_opt_pcg_section, xalmo_opt_trustr_section + penalty_section, xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, & + xalmo_opt_trustr_section, nlmo_opt_trustr_section, nlmo_newton_pcg_section CALL timeset(routineN, handle) @@ -125,12 +125,13 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env) "NLMO_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") + almo_analysis_section => section_vals_get_subs_vals(almo_scf_section, "ANALYSIS") matrix_iterate_section => section_vals_get_subs_vals(almo_scf_section, & "MATRIX_ITERATE") penalty_section => section_vals_get_subs_vals(almo_scf_section, "NLMO_PENALTY") + nlmo_newton_pcg_section => section_vals_get_subs_vals(almo_scf_section, "NLMO_NEWTON_PCG") ! read user input ! common ALMO options @@ -164,6 +165,8 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env) i_val=almo_scf_env%construct_nlmos) CALL section_vals_val_get(almo_scf_section, "VIRTUAL_NLMOS", & l_val=almo_scf_env%virtual_nlmos) + CALL section_vals_val_get(almo_scf_section, "OCCUPIED_NLMOS", & + l_val=almo_scf_env%occupied_nlmos) CALL section_vals_val_get(almo_scf_section, "NLMO_OPERATOR", & i_val=almo_scf_env%nlmo_operator_type) CALL section_vals_val_get(almo_scf_section, "NLMO_COMPACT_FILTER_START", & @@ -318,6 +321,22 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env) "FINAL_DETERMINANT", & r_val=almo_scf_env%opt_nlmo_penalty%final_determinant) + CALL section_vals_val_get(nlmo_newton_pcg_section, "EPS_ERROR", & + r_val=almo_scf_env%opt_nlmo_newton_pcg%eps_error) + CALL section_vals_val_get(nlmo_newton_pcg_section, "MAX_ITER", & + i_val=almo_scf_env%opt_nlmo_newton_pcg%max_iter) + CALL section_vals_val_get(nlmo_newton_pcg_section, "MAX_ITER_OUTER_LOOP", & + i_val=almo_scf_env%opt_nlmo_newton_pcg%max_iter_outer_loop) + CALL section_vals_val_get(nlmo_newton_pcg_section, "START_GRADIENT", & + r_val=almo_scf_env%opt_nlmo_newton_pcg%start_grad) + CALL section_vals_val_get(nlmo_newton_pcg_section, "PRECONDITIONER", & + i_val=almo_scf_env%opt_nlmo_newton_pcg%preconditioner) + + 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", & + i_val=almo_scf_env%almo_analysis%frozen_mo_energy_term) + 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", & @@ -327,10 +346,6 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env) 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(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", & - i_val=almo_scf_env%almo_analysis%frozen_mo_energy_term) !CALL section_vals_val_get(almo_scf_section,"DOMAIN_LAYOUT_AOS",& ! i_val=almo_scf_env%domain_layout_aos) diff --git a/src/almo_scf_methods.F b/src/almo_scf_methods.F index c1b4f7d0d9b..0852e73aa07 100644 --- a/src/almo_scf_methods.F +++ b/src/almo_scf_methods.F @@ -77,10 +77,117 @@ MODULE almo_scf_methods almo_scf_ks_to_ks_xx, & construct_domain_r_down, & xalmo_initial_guess, & + copy_virt2occ, & fill_matrix_with_ones CONTAINS +! ************************************************************************************************** +!> \brief Copy virtual orbitals into occupied orbitals. Skip virtuals that do not fit. +!> \param m_occ ... +!> \param m_virt ... +!> \param first_virtual_of_domain ... +!> \par History +!> 2020.10 created [Rustam Z Khaliullin] +!> \author Rustam Z Khaliullin +! ************************************************************************************************** + SUBROUTINE copy_virt2occ(m_occ, m_virt, first_virtual_of_domain) + + TYPE(dbcsr_type), INTENT(INOUT) :: m_occ + TYPE(dbcsr_type), INTENT(IN) :: m_virt + INTEGER, OPTIONAL, INTENT(IN) :: first_virtual_of_domain + + CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_virt2occ', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, iblock_col, iblock_row, & + iblock_row_size, iblock_col_size, & + nocc_of_block, nvirt_of_block, & + nao_of_block, last_virtual, & + my_first_virtual, & + nvert_blocks2, & + nhori_blocks2, & + nvert_blocks1, & + nhori_blocks1 + LOGICAL :: block_needed + !REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy + REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p, p_new_block + TYPE(dbcsr_iterator_type) :: iter + INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_block_sizes, occ_block_sizes + INTEGER, DIMENSION(:), POINTER :: ao_blk_sizes, occ_blk_sizes + + CALL timeset(routineN, handle) + + my_first_virtual=1 + IF (PRESENT(first_virtual_of_domain)) my_first_virtual=first_virtual_of_domain + CPASSERT(my_first_virtual.GT.0) + + nvert_blocks2 = dbcsr_nblkrows_total(m_occ) + nhori_blocks2 = dbcsr_nblkcols_total(m_occ) + nvert_blocks1 = dbcsr_nblkrows_total(m_virt) + nhori_blocks1 = dbcsr_nblkcols_total(m_virt) + CPASSERT(nvert_blocks1.EQ.nvert_blocks2) + CPASSERT(nhori_blocks1.EQ.nhori_blocks2) + + CALL dbcsr_get_info(m_occ, row_blk_size=ao_blk_sizes) + CALL dbcsr_get_info(m_occ, col_blk_size=occ_blk_sizes) + ALLOCATE (occ_block_sizes(nvert_blocks2), ao_block_sizes(nhori_blocks2)) + occ_block_sizes(:) = occ_blk_sizes(:) + ao_block_sizes(:) = ao_blk_sizes(:) + + CALL dbcsr_set(m_occ, 0.0_dp) + CALL dbcsr_filter(m_occ,1.0_dp) + + CALL dbcsr_work_create(m_occ, work_mutable=.TRUE.) + + CALL dbcsr_iterator_start(iter, m_virt) + 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) + + nocc_of_block = occ_block_sizes(iblock_col) + nao_of_block = ao_block_sizes(iblock_row) + nvirt_of_block = iblock_col_size + CPASSERT(nao_of_block.EQ.iblock_row_size) + + block_needed = .TRUE. + IF ( nocc_of_block*nao_of_block*nvirt_of_block .EQ. 0 ) THEN + block_needed = .FALSE. + ENDIF + IF (my_first_virtual.GT.nvirt_of_block) THEN + block_needed = .FALSE. + ENDIF + + IF (block_needed) THEN + + ! Prepare data + !ALLOCATE (data_copy(nao_of_block, nocc_of_block)) + !data_copy(:, :) = data_p(:, :) + + NULLIFY (p_new_block) + CALL dbcsr_reserve_block2d(m_occ, iblock_row, iblock_col, p_new_block) + CPASSERT(ASSOCIATED(p_new_block)) + p_new_block(:, :) = 0.0_dp + !p_new_block(:, :) = data_p(:, (nocc_of_block + 1):(nocc_of_block + nvirt_of_block)) + last_virtual = MIN(nvirt_of_block, my_first_virtual+nocc_of_block-1) + p_new_block(:, :) = data_p(:, my_first_virtual:last_virtual) + + !DEALLOCATE (data_copy) + + ENDIF + + ENDDO + CALL dbcsr_iterator_stop(iter) + + CALL dbcsr_finalize(m_occ) + + DEALLOCATE (occ_block_sizes) + DEALLOCATE (ao_block_sizes) + + CALL timestop(handle) + + END SUBROUTINE copy_virt2occ + ! ************************************************************************************************** !> \brief Fill all matrix blocks with 1.0_dp !> \param matrix ... diff --git a/src/almo_scf_optimizer.F b/src/almo_scf_optimizer.F index fbfecd80a23..9bf1e0cca8b 100644 --- a/src/almo_scf_optimizer.F +++ b/src/almo_scf_optimizer.F @@ -74,7 +74,7 @@ MODULE almo_scf_optimizer select_row USE input_constants, ONLY: & 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, & + cg_hager_zhang, cg_hestenes_stiefel, cg_liu_storey, cg_polak_ribiere, cg_zero, cg_bfgs_hybrid, & 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_identity USE input_section_types, ONLY: section_vals_get_subs_vals,& @@ -5936,11 +5936,11 @@ END SUBROUTINE compute_preconditioner !> 2015.04 created [Rustam Z Khaliullin] !> \author Rustam Z Khaliullin ! ************************************************************************************************** - SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, & + SUBROUTINE compute_cg_beta(beta, numer, denom, l_index, reset_conjugator, conjugator, & grad, prev_grad, step, prev_step, prev_minus_prec_grad) REAL(KIND=dp), INTENT(INOUT) :: beta - REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: numer, denom + REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: numer, denom, l_index LOGICAL, INTENT(INOUT) :: reset_conjugator INTEGER, INTENT(IN) :: conjugator TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: grad, prev_grad, step, prev_step @@ -5979,7 +5979,8 @@ SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, & IF (PRESENT(numer) .OR. PRESENT(denom)) THEN IF (conjugator .EQ. cg_hestenes_stiefel .OR. & conjugator .EQ. cg_dai_yuan .OR. & - conjugator .EQ. cg_hager_zhang) THEN + conjugator .EQ. cg_hager_zhang .OR. & + conjugator .EQ. cg_bfgs_hybrid) THEN CPABORT("cannot return numer/denom") ENDIF ENDIF @@ -6034,6 +6035,16 @@ SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, & CALL dbcsr_dot(prev_step(i), grad(i), num3) my_numer2 = my_numer2 + num2 my_numer3 = my_numer3 + num3 + CASE (cg_bfgs_hybrid) + CALL dbcsr_copy(m_tmp_no_1, grad(i)) + CALL dbcsr_add(m_tmp_no_1, prev_grad(i), & + 1.0_dp, -1.0_dp) + CALL dbcsr_dot(m_tmp_no_1, grad(i), num) + CALL dbcsr_dot(prev_grad(i), prev_step(i), den) + CALL dbcsr_dot(grad(i), grad(i), num2) + CALL dbcsr_dot(grad(i), prev_step(i), num3) + my_numer2 = my_numer2 + num2 + my_numer3 = my_numer3 + num3 CASE (cg_zero) num = 0.0_dp den = 1.0_dp @@ -6052,12 +6063,16 @@ SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, & SELECT CASE (conjugator) CASE (cg_hestenes_stiefel, cg_dai_yuan) beta = -1.0_dp*my_numer/my_denom - CASE (cg_fletcher_reeves, cg_polak_ribiere, cg_fletcher, cg_liu_storey) + CASE (cg_fletcher_reeves, cg_polak_ribiere, & + cg_fletcher, cg_liu_storey) beta = my_numer/my_denom CASE (cg_hager_zhang) kappa = -2.0_dp*my_numer/my_denom tau = -1.0_dp*my_numer2/my_denom beta = tau - kappa*my_numer3/my_denom + CASE (cg_bfgs_hybrid) + beta = MAX(0.0_dp, MIN(-my_numer/my_denom, -my_numer2/my_denom)) + l_index = 1.0_dp + beta*(my_numer3/my_numer2) CASE (cg_zero) beta = 0.0_dp CASE DEFAULT @@ -9591,7 +9606,7 @@ SUBROUTINE trust_r_report(unit_nr, iter_type, iteration, radius, & CASE (1) - WRITE (unit_nr, '(T2,A6,A5,I6,F22.10,ES10.2,T67,ES7.0,F6.1)') & + WRITE (unit_nr, '(T2,A6,A5,I6,F22.10,ES14.5,T67,ES7.0,F10.5)') & iter_type_str, & iter_status, & iteration, & @@ -9602,11 +9617,12 @@ SUBROUTINE trust_r_report(unit_nr, iter_type, iteration, radius, & CASE (2) - WRITE (unit_nr, '(T2,A6,A5,I6,F22.10,ES10.2,ES10.2,F6.1,ES7.0,F6.1)') & + WRITE (unit_nr, '(T2,A6,A5,I6,F22.10,ES14.5,ES10.2,ES10.2,F6.1,ES5.0,F10.5)') & iter_type_str, & iter_status, & iteration, & loss, & + grad_norm, & delta_loss, predicted_reduction, rho, & ! distinct radius, & time diff --git a/src/almo_scf_types.F b/src/almo_scf_types.F index 94adc09b0be..3d44d525849 100644 --- a/src/almo_scf_types.F +++ b/src/almo_scf_types.F @@ -18,7 +18,7 @@ MODULE almo_scf_types USE domain_submatrix_types, ONLY: domain_map_type,& domain_submatrix_type USE input_constants, ONLY: & - cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, & + cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, cg_bfgs_hybrid, & cg_liu_storey, cg_polak_ribiere, cg_zero, optimizer_diis, optimizer_pcg, optimizer_trustr, & trustr_cauchy, trustr_dogleg, trustr_steihaug, xalmo_prec_domain, xalmo_prec_full, & xalmo_prec_identity @@ -51,6 +51,16 @@ MODULE almo_scf_types END TYPE penalty_type + TYPE nlmo_newton_pcg_type + + REAL(KIND=dp) :: eps_error, start_grad + + INTEGER :: preconditioner, & ! preconditioner type + max_iter, & + max_iter_outer_loop + + END TYPE nlmo_newton_pcg_type + ! almo-based electronic structure analysis TYPE almo_analysis_type @@ -212,9 +222,11 @@ MODULE almo_scf_types ! NLMO-related INTEGER :: construct_nlmos LOGICAL :: virtual_nlmos + LOGICAL :: occupied_nlmos INTEGER :: nlmo_operator_type REAL(KIND=dp) :: nlmo_compactification_filter_start TYPE(penalty_type) :: opt_nlmo_penalty + TYPE(nlmo_newton_pcg_type) :: opt_nlmo_newton_pcg TYPE(optimizer_options_type) :: opt_nlmo_trustr TYPE(optimizer_options_type) :: opt_nlmo_pcg @@ -375,6 +387,8 @@ MODULE almo_scf_types TYPE(optimizer_options_type) :: opt_xalmo_trustr TYPE(optimizer_options_type) :: opt_xalmo_newton_pcg_solver TYPE(optimizer_options_type) :: opt_k_pcg + TYPE(optimizer_options_type) :: opt_nlmo_pcg_newton_solver + TYPE(optimizer_options_type) :: opt_nlmo_trustr_newton_solver ! keywords that control electron delocalization treatment ! RZK-warning: many of these varibles should be collected @@ -497,6 +511,8 @@ SUBROUTINE print_optimizer_options(optimizer, unit_nr) conj_string = "Dai-Yuan" CASE (cg_hager_zhang) conj_string = "Hager-Zhang" + CASE (cg_bfgs_hybrid) + conj_string = "BFGS-Hybrid" END SELECT WRITE (unit_nr, '(T4,A,T48,A33)') "conjugator:", TRIM(conj_string) diff --git a/src/input_constants.F b/src/input_constants.F index 39f18aa4a90..970f2af66d0 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -936,7 +936,8 @@ MODULE input_constants cg_fletcher = 4, & cg_liu_storey = 5, & cg_dai_yuan = 6, & - cg_hager_zhang = 7 + cg_hager_zhang = 7, & + cg_bfgs_hybrid = 8 INTEGER, PARAMETER, PUBLIC :: trustr_steihaug = 1, & trustr_cauchy = 2, & trustr_dogleg = 3 diff --git a/src/input_cp2k_almo.F b/src/input_cp2k_almo.F index a0f843519e9..4f13a14f14e 100644 --- a/src/input_cp2k_almo.F +++ b/src/input_cp2k_almo.F @@ -19,8 +19,8 @@ MODULE input_cp2k_almo almo_deloc_xalmo_1diag, almo_deloc_xalmo_scf, almo_deloc_xalmo_x, almo_frz_crystal, & 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, & + cg_liu_storey, cg_polak_ribiere, cg_zero, cg_bfgs_hybrid, 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_identity, xalmo_prec_lbfgs, xalmo_trial_r0_out, xalmo_trial_simplex @@ -231,6 +231,13 @@ SUBROUTINE create_almo_scf_section(section) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="OCCUPIED_NLMOS", & + description="Decide if localizing occupied orbitals", & + usage="OCCUPIED_NLMOS .TRUE.", default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="NLMO_COMPACT_FILTER_START", & description="Set orbital coefficients with absolute value smaller than this value to zero.", & usage="NLMO_COMPACT_FILTER_START 1.e-6", default_r_val=-1.0_dp) @@ -583,6 +590,11 @@ SUBROUTINE create_almo_scf_section(section) CALL section_add_subsection(section, subsection) CALL section_release(subsection) + NULLIFY (subsection) + CALL create_nlmo_newton_pcg_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + NULLIFY (subsection) CALL create_matrix_iterate_section(subsection) CALL section_add_subsection(section, subsection) @@ -636,6 +648,7 @@ RECURSIVE SUBROUTINE create_optimizer_section(section, optimizer_id) 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.) + NULLIFY (subsection) optimizer_type = optimizer_pcg CASE (optimizer_xalmo_pcg) CALL section_create(section, __LOCATION__, name="XALMO_OPTIMIZER_PCG", & @@ -654,6 +667,7 @@ RECURSIVE SUBROUTINE create_optimizer_section(section, optimizer_id) "approaches. An iterative conjugate-gradient approach is "// & "used and controled by the inner loop", & n_keywords=10, n_subsections=0, repeats=.FALSE.) + NULLIFY (subsection) optimizer_type = optimizer_trustr CASE (optimizer_xalmo_trustr) CALL section_create(section, __LOCATION__, name="XALMO_OPTIMIZER_TRUSTR", & @@ -780,14 +794,15 @@ RECURSIVE SUBROUTINE create_optimizer_section(section, optimizer_id) usage="CONJUGATOR POLAK_RIBIERE", & default_i_val=cg_hager_zhang, & enum_c_vals=s2a("ZERO", "POLAK_RIBIERE", "FLETCHER_REEVES", & - "HESTENES_STIEFEL", "FLETCHER", "LIU_STOREY", "DAI_YUAN", "HAGER_ZHANG"), & + "HESTENES_STIEFEL", "FLETCHER", "LIU_STOREY", "DAI_YUAN", "HAGER_ZHANG", & + "BFGS_HYBRID"), & enum_desc=s2a("Steepest descent", "Polak and Ribiere", & "Fletcher and Reeves", "Hestenes and Stiefel", & "Fletcher (Conjugate descent)", "Liu and Storey", & - "Dai and Yuan", "Hager and Zhang"), & + "Dai and Yuan", "Hager and Zhang", "CG-BFGS Hybrid"), & enum_i_vals=(/cg_zero, cg_polak_ribiere, cg_fletcher_reeves, & cg_hestenes_stiefel, cg_fletcher, cg_liu_storey, & - cg_dai_yuan, cg_hager_zhang/)) + cg_dai_yuan, cg_hager_zhang, cg_bfgs_hybrid/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -940,6 +955,73 @@ SUBROUTINE create_penalty_section(section) END SUBROUTINE create_penalty_section +! ************************************************************************************************** +!> \brief The section controls an iterative solver of the Newton-Raphson linear equation. +!> \param section ... +!> \par History +!> 2020.07 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE create_nlmo_newton_pcg_section(section) + + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + NULLIFY (section) + + CALL section_create(section, __LOCATION__, name="NLMO_NEWTON_PCG", & + description="Controls an iterative solver of the Newton-Raphson linear equation.", & + n_keywords=5, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER", & + description="Maximum number of iterations", & + usage="MAX_ITER 100", default_i_val=20) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EPS_ERROR", & + description="Target value of the MAX norm of the error", & + usage="EPS_ERROR 1.E-6", default_r_val=1.0E-5_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER_OUTER_LOOP", & + description="Maximum number of iterations in the outer loop. "// & + "Use the outer loop to update the preconditioner and reset the conjugator. "// & + "This can speed up convergence significantly.", & + usage="MAX_ITER 10", default_i_val=0) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="START_GRADIENT", & + description="Avoid initial possible negative Hessian stage by"// & + "swiching off the newton solver inside an iteration", & + usage="START_GRADIENT 10.0", default_r_val=0.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="PRECONDITIONER", & + description="Select a preconditioner for the conjugate gradient optimization", & + usage="PRECONDITIONER DOMAIN", & + default_i_val=xalmo_prec_domain, & + enum_c_vals=s2a("NONE", "DEFAULT", "DOMAIN", "FULL", "LBFGS"), & + enum_desc=s2a("Do not use preconditioner", & + "Same as DOMAIN preconditioner", & + "Invert preconditioner domain-by-domain."// & + " The main component of the linear scaling algorithm", & + "Solve linear equations step=-H.grad on the entire space", & + "Use limited-memory BFGS method"), & + enum_i_vals=(/xalmo_prec_identity, xalmo_prec_domain, & + xalmo_prec_domain, xalmo_prec_full, xalmo_prec_lbfgs/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_nlmo_newton_pcg_section + ! ************************************************************************************************** !> \brief The section controls electronic structure analysis based on ALMOs !> \param section ... diff --git a/src/nlmo_methods.F b/src/nlmo_methods.F index 28fa3c39fb1..27d61ac8761 100644 --- a/src/nlmo_methods.F +++ b/src/nlmo_methods.F @@ -16,20 +16,38 @@ MODULE nlmo_methods lbfgs_seed USE almo_scf_methods, ONLY: fill_matrix_with_ones USE almo_scf_qs, ONLY: matrix_qs_to_almo - USE almo_scf_types, ONLY: optimizer_options_type + USE almo_scf_optimizer, ONLY: compute_cg_beta + USE almo_scf_types, ONLY: optimizer_options_type, & + almo_scf_env_type USE cell_types, ONLY: cell_type - USE cp_log_handling, ONLY: cp_to_string + USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& + cp_dbcsr_cholesky_invert,& + cp_dbcsr_cholesky_restore + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_unit_nr,& + cp_logger_type,& + cp_to_string + USE cp_para_types, ONLY: cp_para_env_type + USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& + dbcsr_deallocate_matrix_set USE dbcsr_api, ONLY: & - dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_get_diag, & - dbcsr_get_info, dbcsr_hadamard_product, dbcsr_multiply, dbcsr_p_type, dbcsr_release, & - dbcsr_scale, dbcsr_set, dbcsr_set_diag, dbcsr_type, dbcsr_type_no_symmetry - USE input_constants, ONLY: cg_zero,& - op_loc_berry,& - op_loc_pipek,& - xalmo_prec_dbfgs,& - xalmo_prec_full,& - xalmo_prec_identity,& - xalmo_prec_lbfgs + dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, & + dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_finalize, & + dbcsr_frobenius_norm, dbcsr_func_dtanh, dbcsr_func_inverse, dbcsr_func_tanh, & + dbcsr_function_of_elements, dbcsr_get_block_p, dbcsr_get_diag, dbcsr_get_info, & + dbcsr_hadamard_product, 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_norm, dbcsr_norm_maxabsnorm, & + dbcsr_p_type, dbcsr_print_block_sum, dbcsr_release, dbcsr_reserve_block2d, dbcsr_scale, & + dbcsr_set, dbcsr_set_diag, dbcsr_triu, dbcsr_type, dbcsr_type_no_symmetry, & + dbcsr_work_create, dbcsr_print + USE input_constants, ONLY: & + 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, cg_bfgs_hybrid, & + 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_identity, op_loc_berry,& + op_loc_pipek, xalmo_prec_dbfgs, xalmo_prec_lbfgs USE iterate_matrix, ONLY: determinant,& invert_Hotelling USE kinds, ONLY: dp @@ -66,7 +84,9 @@ MODULE nlmo_methods nlmo_env_minus_hessian_inv_apply, & nlmo_env_invert_hessian, & nlmo_env_hessian_apply, & - nlmo_env_predicted_reduction + nlmo_newton_grad_to_step, & + nlmo_env_predicted_reduction, & + trust_step_subproblem CONTAINS @@ -127,8 +147,12 @@ SUBROUTINE nlmo_env_init(nlmo_env, qs_env, loc_operator, m_templateNN, & nlmo_env%weights = 0.0_dp CALL initialize_weights(cell, nlmo_env%weights) - ALLOCATE (op_sm_set_qs(2, dim_op)) - ALLOCATE (nlmo_env%op_sm_set(2, dim_op)) + NULLIFY (op_sm_set_qs) + CALL dbcsr_allocate_matrix_set(op_sm_set_qs, 2, dim_op) + NULLIFY (nlmo_env%op_sm_set) + CALL dbcsr_allocate_matrix_set(nlmo_env%op_sm_set, 2, dim_op) + !ALLOCATE (op_sm_set_qs(2, dim_op)) + !ALLOCATE (nlmo_env%op_sm_set(2, dim_op)) ! allocate and initialize all matrices DO idim0 = 1, dim_op ! this loop is over miller ind @@ -138,6 +162,8 @@ SUBROUTINE nlmo_env_init(nlmo_env, qs_env, loc_operator, m_templateNN, & nlmo_env%op_sm_set(reim, idim0)%matrix) ALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + CALL dbcsr_create(op_sm_set_qs(reim, idim0)%matrix, & + template=qs_matrix_s(1)%matrix) CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, & qs_matrix_s(1)%matrix, & name="QS_NLMO_"// & @@ -146,6 +172,8 @@ SUBROUTINE nlmo_env_init(nlmo_env, qs_env, loc_operator, m_templateNN, & CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp) ALLOCATE (nlmo_env%op_sm_set(reim, idim0)%matrix) + CALL dbcsr_create(nlmo_env%op_sm_set(reim, idim0)%matrix, & + template=m_templateNN) CALL dbcsr_copy(nlmo_env%op_sm_set(reim, idim0)%matrix, & m_templateNN, & name="NLMO_"// & @@ -166,10 +194,11 @@ SUBROUTINE nlmo_env_init(nlmo_env, qs_env, loc_operator, m_templateNN, & CALL matrix_qs_to_almo(op_sm_set_qs(reim, idim0)%matrix, & nlmo_env%op_sm_set(reim, idim0)%matrix, & distr_type_AOs, .FALSE.) - DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix) + !DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix) ENDDO ENDDO - DEALLOCATE (op_sm_set_qs) + !DEALLOCATE (op_sm_set_qs) + CALL dbcsr_deallocate_matrix_set(op_sm_set_qs) CASE (op_loc_pipek) @@ -355,17 +384,19 @@ SUBROUTINE nlmo_env_release(nlmo_env) CALL dbcsr_release(nlmo_env%m_B0(reim, idim0, ispin)) ENDDO - SELECT CASE (nlmo_env%loc_operator) - CASE (op_loc_berry) - DEALLOCATE (nlmo_env%op_sm_set(reim, idim0)%matrix) - END SELECT + !SELECT CASE (nlmo_env%loc_operator) + !CASE (op_loc_berry) + ! DEALLOCATE (nlmo_env%op_sm_set(reim, idim0)%matrix) + !END SELECT ENDDO ENDDO SELECT CASE (nlmo_env%loc_operator) CASE (op_loc_berry) - DEALLOCATE (nlmo_env%op_sm_set) + !DEALLOCATE (nlmo_env%op_sm_set) + CALL dbcsr_deallocate_matrix_set(nlmo_env%op_sm_set) + END SELECT DEALLOCATE (nlmo_env%m_B0) @@ -399,8 +430,10 @@ SUBROUTINE nlmo_env_set_flags(nlmo_env, optimizer) ELSE IF (nlmo_env%hessian_type .EQ. xalmo_prec_dbfgs) THEN nlmo_env%d_bfgs = .TRUE. ENDIF - IF (nlmo_env%l_bfgs .AND. (optimizer%conjugator .NE. cg_zero)) THEN - CPABORT("Cannot use conjugators with BFGS") + + IF (nlmo_env%l_bfgs .AND. (optimizer%conjugator .NE. cg_zero .AND. & + optimizer%conjugator .NE. cg_bfgs_hybrid)) THEN + CPABORT("Cannot use conjugators with BFGS other than hybrid") ENDIF END SUBROUTINE nlmo_env_set_flags @@ -591,6 +624,24 @@ SUBROUTINE nlmo_env_mainvar_to_aux(nlmo_env) INTEGER :: ispin, nocc, nspins, para_group REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diagonal TYPE(dbcsr_type) :: tm_OO1 + REAL(KIND=dp) :: t_invert_Hotelling, & + t_dbcsr_multi + INTEGER :: unit_nr, handle + TYPE(cp_logger_type), POINTER :: logger + CHARACTER(len=*), PARAMETER :: routineN = 'nlmo_env_mainvar_to_aux', & + routineP = moduleN//':'//routineN + + 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 + + t_invert_Hotelling = 0.0_dp + t_dbcsr_multi = 0.0_dp nspins = SIZE(nlmo_env%m_theta) @@ -612,6 +663,7 @@ SUBROUTINE nlmo_env_mainvar_to_aux(nlmo_env) 0.0_dp, & tm_OO1, & filter_eps=nlmo_env%eps_filter) + CALL dbcsr_set(nlmo_env%m_sig_sqrti_ii(ispin), 0.0_dp) CALL dbcsr_add_on_diag(nlmo_env%m_sig_sqrti_ii(ispin), 1.0_dp) CALL dbcsr_multiply("T", "N", 1.0_dp, & @@ -673,6 +725,8 @@ SUBROUTINE nlmo_env_mainvar_to_aux(nlmo_env) ENDDO ! ispin + CALL timestop(handle) + END SUBROUTINE nlmo_env_mainvar_to_aux ! ************************************************************************************************** @@ -699,13 +753,24 @@ SUBROUTINE nlmo_env_loss_function(nlmo_env, loss_components, overlap_determinant penal_function_ispin REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: reim_diag, z2 TYPE(dbcsr_type) :: tempOccOcc1, tempOccOcc2 + INTEGER :: unit_nr + TYPE(cp_logger_type), POINTER :: logger 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(nlmo_env%m_theta_normalized) nlmo_env%localization_component = 0.0_dp nlmo_env%orthogonalization_component = 0.0_dp + penal_function_ispin = 0.0_dp + det1 = 0.0_dp DO ispin = 1, nspins CALL dbcsr_get_info(nlmo_env%m_theta_normalized(ispin), & @@ -743,7 +808,7 @@ SUBROUTINE nlmo_env_loss_function(nlmo_env, loss_components, overlap_determinant 0.0_dp, tempOccOcc2, & retain_sparsity=.TRUE.) - reim_diag = 0.0_dp + 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(:) @@ -760,6 +825,10 @@ SUBROUTINE nlmo_env_loss_function(nlmo_env, loss_components, overlap_determinant fval = nlmo_env%weights(idim0) - nlmo_env%weights(idim0)*SQRT(ABS(z2(ielem))) END SELECT local_function_ispin = local_function_ispin + fval + !IF (unit_nr > 0) THEN + ! WRITE (unit_nr, '(T2,3A19,I6,I6,F23.10)') & + ! "Atom", "No. Occ", "Localization:", idim0, ielem, fval + !ENDIF ENDDO ENDDO ! end loop over idim0 @@ -812,10 +881,25 @@ SUBROUTINE nlmo_env_loss_gradient(nlmo_env) 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 + REAL(KIND=dp) :: t_grad_loc, t_grad_pen, t_grad_norm, & + t_dbcsr_multi_loc + INTEGER :: unit_nr + TYPE(cp_logger_type), POINTER :: logger 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(nlmo_env%m_theta_normalized) + t_grad_loc = 0.0_dp + t_grad_pen = 0.0_dp + t_grad_norm = 0.0_dp + t_dbcsr_multi_loc = 0.0_dp DO ispin = 1, nspins @@ -835,6 +919,7 @@ SUBROUTINE nlmo_env_loss_gradient(nlmo_env) CALL dbcsr_get_info(nlmo_env%m_siginv(ispin), 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 @@ -896,7 +981,6 @@ SUBROUTINE nlmo_env_loss_gradient(nlmo_env) ENDDO ! end loop over idim0 DEALLOCATE (z2) - ! add gradient of the penalty functional log[det(sigma)] ! G = 2*prefactor*sigma0.a_norm.sigma_inv !CALL dbcsr_norm(nlmo_env%m_sigma0_thetanorm_siginv(ispin), & @@ -1156,15 +1240,15 @@ END SUBROUTINE nlmo_env_invert_hessian !> \param m_in ... !> \param m_out ... !> \par History -!> 2020.02 created [Rustam Z Khaliullin] -!> \author Rustam Z Khaliullin +!> 2020.02 created [Rustam Z Khaliullin, Ziling Luo] +!> \author Rustam Z Khaliullin, Ziling Luo ! ************************************************************************************************** SUBROUTINE nlmo_env_hessian_apply(nlmo_env, m_in, m_out) - TYPE(nlmo_env_type), INTENT(INOUT) :: nlmo_env + TYPE(nlmo_env_type), INTENT(IN) :: nlmo_env + TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_out TYPE(dbcsr_type), DIMENSION(:), INTENT(IN), & OPTIONAL :: m_in - TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_out INTEGER :: dim0, idim0, ispin, nspins, reim LOGICAL :: d_bfgs, l_bfgs @@ -1177,16 +1261,17 @@ SUBROUTINE nlmo_env_hessian_apply(nlmo_env, m_in, m_out) d_bfgs = nlmo_env%d_bfgs l_bfgs = nlmo_env%l_bfgs - IF (d_bfgs .OR. l_bfgs) THEN + IF (d_bfgs) THEN ! RZK-hessian: unclear what needs to be done here ! Do we use BFGS for the PCG/TRUSTR. Within TRUSTR: is it a subproblem solver? - CPABORT("This Hessian is NYI for BFGS (sub)problem solver") + !CPABORT("This Hessian is NYI for BFGS (sub)problem solver") ELSE ! non-BFGS - SELECT CASE (nlmo_env%hessian_type) - CASE (xalmo_prec_full) + !SELECT CASE (nlmo_env%hessian_type) + !CASE (xalmo_prec_full) + IF (nlmo_env%hessian_type .EQ. xalmo_prec_full .OR. l_bfgs) THEN DO ispin = 1, nspins @@ -1247,7 +1332,7 @@ SUBROUTINE nlmo_env_hessian_apply(nlmo_env, m_in, m_out) 1.0_dp, tm_oo_1, & filter_eps=nlmo_env%eps_filter) - ! Accumulate: (-1)*T.(tr(m_in).T)_diag.N, where T=sigma0.A + ! Accumulate: (-2)*T.(tr(m_in).T)_diag.N, where T=sigma0.A CALL dbcsr_multiply("T", "N", 1.0_dp, & m_in(ispin), & nlmo_env%m_sigma0_thetanorm(ispin), & @@ -1262,7 +1347,7 @@ SUBROUTINE nlmo_env_hessian_apply(nlmo_env, m_in, m_out) tm_oo_2, & 0.0_dp, tm_oo_3, & filter_eps=nlmo_env%eps_filter) - CALL dbcsr_multiply("N", "N", -1.0_dp, & + CALL dbcsr_multiply("N", "N", -2.0_dp, & tm_oo_3, & nlmo_env%m_sig_sqrti_ii(ispin), & 1.0_dp, tm_oo_1, & @@ -1453,14 +1538,16 @@ SUBROUTINE nlmo_env_hessian_apply(nlmo_env, m_in, m_out) ENDDO - CASE (xalmo_prec_identity) + !CASE (xalmo_prec_identity) + ELSE IF (nlmo_env%hessian_type .EQ. xalmo_prec_identity) THEN ! Hessian is set to identity by user DO ispin = 1, nspins CALL dbcsr_copy(m_out(ispin), m_in(ispin)) ENDDO ! ispin - END SELECT + !END SELECT + ENDIF ENDIF ! hessian type @@ -1553,39 +1640,71 @@ END SUBROUTINE nlmo_env_minus_hessian_inv_apply ! ************************************************************************************************** !> \brief This subroutine is obsolete. Do not call !> \param nlmo_env ... +!> \param almo_scf_env ... !> \param grad ... !> \param step ... !> \param prev_grad ... !> \param prev_m_theta ... !> \param iteration ... +!> \param identity_on ... !> \par History !> 2020.02 created [Rustam Z Khaliullin] !> \author Rustam Z Khaliullin ! ************************************************************************************************** - SUBROUTINE nlmo_env_loss_curvature_grad_to_step(nlmo_env, grad, step, & + SUBROUTINE nlmo_env_loss_curvature_grad_to_step(nlmo_env, almo_scf_env, grad, step, & prev_grad, prev_m_theta, & - iteration) + iteration, identity_on) !RZK-critical: refactor using other hessian methods TYPE(nlmo_env_type), INTENT(INOUT) :: nlmo_env - TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: grad + TYPE(almo_scf_env_type), INTENT(IN) :: almo_scf_env TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: step - TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: prev_grad, prev_m_theta - INTEGER :: iteration + TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: grad, prev_grad, & + prev_m_theta + INTEGER, INTENT(IN) :: iteration + LOGICAL, INTENT(IN) :: identity_on CHARACTER(len=*), PARAMETER :: routineN = 'nlmo_env_loss_curvature_grad_to_step', & routineP = moduleN//':'//routineN - INTEGER :: ispin, nspins + INTEGER :: ispin, nspins, unit_nr LOGICAL :: d_bfgs, l_bfgs - REAL(KIND=dp) :: bfgs_rho, bfgs_sum - TYPE(dbcsr_type) :: bfgs_s, bfgs_y, tempOccOcc1, & - tempOccOcc2, tempOccOcc3 + REAL(KIND=dp) :: bfgs_rho, bfgs_sum, check_norm + TYPE(dbcsr_type) :: tempOccOcc1, & + tempOccOcc2, & + tempOccOcc3, & + bfgs_s, & + bfgs_y + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: m_Hstep + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: check_max_norm + TYPE(cp_logger_type), POINTER :: logger + + ! get a useful 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(nlmo_env%m_theta) d_bfgs = nlmo_env%d_bfgs l_bfgs = nlmo_env%l_bfgs + + ALLOCATE (m_Hstep(nspins)) + + ALLOCATE (check_max_norm(nspins)) + + DO ispin = 1, nspins + + ! init matrices + CALL dbcsr_create(m_Hstep(ispin), & + template=nlmo_env%m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + + ENDDO !ispin + ! if available use second derivative info - bfgs, hessian, preconditioner IF (nlmo_env%hessian_type .EQ. xalmo_prec_identity) THEN ! no second derivatives @@ -1597,6 +1716,39 @@ SUBROUTINE nlmo_env_loss_curvature_grad_to_step(nlmo_env, grad, step, & ENDDO ! ispin + ELSE IF (nlmo_env%hessian_type .EQ. xalmo_prec_full) THEN + + IF (identity_on) THEN + + DO ispin = 1, nspins + + CALL dbcsr_copy(step(ispin), grad(ispin)) + CALL dbcsr_scale(step(ispin), -1.0_dp) + + ENDDO ! ispin + + ELSE + + CALL nlmo_newton_grad_to_step(nlmo_env, almo_scf_env, step) + ! double check if newton solver actually works + ! calculate Norm(H*Step + G) + CALL nlmo_env_hessian_apply(nlmo_env, step, m_Hstep) + DO ispin = 1, nspins + + CALL dbcsr_add(m_Hstep(ispin), grad(ispin), 1.0_dp, 1.0_dp) + CALL dbcsr_norm(m_Hstep(ispin), dbcsr_norm_maxabsnorm, & + norm_scalar=check_max_norm(ispin)) + + ENDDO ! ispin + check_norm = MAXVAL(check_max_norm) + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T6,A55, F15.10)') & + "check if the newton solver correctly predict the step:", & + check_norm + ENDIF + + ENDIF + ELSE ! use second derivatives ! compute and invert hessian/precond? @@ -1722,8 +1874,234 @@ SUBROUTINE nlmo_env_loss_curvature_grad_to_step(nlmo_env, grad, step, & ENDIF ! second derivative type fork + DO ispin = 1, nspins + CALL dbcsr_release(m_Hstep(ispin)) + ENDDO ! ispin + + DEALLOCATE (m_Hstep) + + END SUBROUTINE nlmo_env_loss_curvature_grad_to_step +! ***************************************************************************** +!> \brief computes the step matrix from the gradient and Hessian using +!> the Newton-Raphson method +!> \param nlmo_env ... +!> \param almo_scf_env ... +!> \param m_delta ... +!> \par History +!> 2020.06 created [Ziling Luo] +!> \author Ziling Luo +! ************************************************************************************************** + SUBROUTINE nlmo_newton_grad_to_step(nlmo_env, almo_scf_env, m_delta) + + + TYPE(nlmo_env_type), INTENT(IN) :: nlmo_env + TYPE(almo_scf_env_type), INTENT(IN) :: almo_scf_env + TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: m_delta + + CHARACTER(len=*), PARAMETER :: routineN = 'newton_grad_to_step', & + routineP = moduleN//':'//routineN + CHARACTER(LEN=20) :: iter_type + INTEGER :: handle, unit_nr, ispin, iteration, & + nspins, outer_iteration + LOGICAL :: converged, outer_prepare_to_exit, & + prepare_to_exit, reset_conjugator, & + use_preconditioner + REAL(KIND=dp) :: alpha, beta, denom, denom_ispin, & + numer, numer_ispin, residue_norm, & + t1, t2, obj_func + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: residue_max_norm + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: m_Hstep, m_step, & + m_residue, m_residue_prev + + + CALL timeset(routineN, handle) + + ! get a useful 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(nlmo_env%m_theta) + reset_conjugator = .TRUE. + use_preconditioner = (almo_scf_env%opt_nlmo_newton_pcg%preconditioner .NE. xalmo_prec_identity) + + ! allocate matrices + ALLOCATE (m_residue(nspins)) + ALLOCATE (m_residue_prev(nspins)) + ALLOCATE (m_step(nspins)) + ALLOCATE (m_Hstep(nspins)) + + ALLOCATE (residue_max_norm(nspins)) + + ! initiate objects before iterations + DO ispin = 1, nspins + + ! init matrices + CALL dbcsr_create(m_residue(ispin), & + template=nlmo_env%m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_residue_prev(ispin), & + template=nlmo_env%m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_step(ispin), & + template=nlmo_env%m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(m_Hstep(ispin), & + template=nlmo_env%m_theta(ispin), & + matrix_type=dbcsr_type_no_symmetry) + + ! initial guess + CALL dbcsr_set(m_delta(ispin), 0.0_dp) + ! r_0 = Ax-b = g + CALL dbcsr_copy(m_residue(ispin), nlmo_env%grad(ispin)) + ! p_0 = -r_0 = -g + CALL dbcsr_copy(m_step(ispin), nlmo_env%grad(ispin)) + CALL dbcsr_scale(m_step(ispin), -1.0_dp) + + IF (use_preconditioner) THEN + + CPABORT("Preconditioner is currently not implemented") + + ENDIf ! do not use preconditioner + + ENDDO !ispin + + ! start the outer SCF loop + outer_prepare_to_exit = .FALSE. + outer_iteration = 0 + residue_norm = 0.0_dp + + DO + + ! start the inner SCF loop + prepare_to_exit = .FALSE. + converged = .FALSE. + iteration = 0 + + DO + + CALL nlmo_env_hessian_apply(nlmo_env, m_delta, m_Hstep) + + obj_func = 0.0_dp + ! Compute current localization function value + DO ispin = 1, nspins + + CALL dbcsr_dot(m_delta(ispin), m_Hstep(ispin), numer_ispin) + CALL dbcsr_dot(m_delta(ispin), nlmo_env%grad(ispin), denom_ispin) + + obj_func = obj_func + 0.5_dp * numer_ispin + denom_ispin + + ENDDO !ispin + + CALL nlmo_env_hessian_apply(nlmo_env, m_step, m_Hstep) + + ! alpha is computed outside the spin loop + numer = 0.0_dp + denom = 0.0_dp + DO ispin = 1, nspins + + CALL dbcsr_dot(m_residue(ispin), m_residue(ispin), numer_ispin) + CALL dbcsr_dot(m_step(ispin), m_Hstep(ispin), denom_ispin) + + numer = numer + numer_ispin + denom = denom + denom_ispin + + ENDDO !ispin + + alpha = numer/denom + + DO ispin = 1, nspins + + ! update the variable + CALL dbcsr_add(m_delta(ispin), m_step(ispin), 1.0_dp, alpha) + CALL dbcsr_copy(m_residue_prev(ispin), m_residue(ispin)) + CALL dbcsr_add(m_residue(ispin), m_Hstep(ispin), & + 1.0_dp, 1.0_dp*alpha) + CALL dbcsr_norm(m_residue(ispin), dbcsr_norm_maxabsnorm, & + norm_scalar=residue_max_norm(ispin)) + + ENDDO ! ispin + + ! check convergence and other exit criteria + residue_norm = MAXVAL(residue_max_norm) + converged = (residue_norm .LT. almo_scf_env%opt_nlmo_newton_pcg%eps_error) + IF (converged .OR. (iteration .GE. almo_scf_env%opt_nlmo_newton_pcg%max_iter)) THEN + prepare_to_exit = .TRUE. + ENDIF + + IF (.NOT. prepare_to_exit) THEN + + ! compute the conjugation coefficient - beta + CALL compute_cg_beta( & + beta=beta, & + reset_conjugator=reset_conjugator, & + conjugator=cg_fletcher, & + grad=m_residue, & + prev_grad=m_residue_prev, & + step=m_residue, & + prev_step=m_residue_prev) + + DO ispin = 1, nspins + + ! conjugate the step direction + CALL dbcsr_add(m_step(ispin), m_residue(ispin), beta, -1.0_dp) + + ENDDO !ispin + + ENDIF ! not.prepare_to_exit + + IF (unit_nr > 0) THEN + iter_type = TRIM("NR STEP") + WRITE (unit_nr, '(T6,A9,I6,F14.5,F14.5,F15.10,F9.2)') & + iter_type, iteration, & + alpha, obj_func, residue_norm, & + t2 - t1 + ENDIF + + iteration = iteration + 1 + IF (prepare_to_exit) EXIT + + ENDDO ! inner loop + + IF (converged .OR. (outer_iteration .GE. almo_scf_env%opt_nlmo_newton_pcg%max_iter_outer_loop)) THEN + outer_prepare_to_exit = .TRUE. + ENDIF + + outer_iteration = outer_iteration + 1 + IF (outer_prepare_to_exit) EXIT + + ENDDO ! outer loop + + + ! clean up + DO ispin = 1, nspins + CALL dbcsr_release(m_residue(ispin)) + CALL dbcsr_release(m_residue_prev(ispin)) + CALL dbcsr_release(m_step(ispin)) + CALL dbcsr_release(m_Hstep(ispin)) + ENDDO !ispin + DEALLOCATE (m_residue) + DEALLOCATE (m_residue_prev) + DEALLOCATE (m_step) + DEALLOCATE (residue_max_norm) + DEALLOCATE (m_Hstep) + + IF (.NOT. converged) THEN + CPABORT("Optimization not converged!") + ENDIF + + ! check that the step satisfies H.step=-grad + + CALL timestop(handle) + + END SUBROUTINE nlmo_newton_grad_to_step + ! ************************************************************************************************** !> \brief ... !> \param nlmo_env ... @@ -1778,6 +2156,42 @@ SUBROUTINE nlmo_env_deallocate_all(nlmo_env) END SUBROUTINE nlmo_env_deallocate_all +! ************************************************************************************************** +!> \brief Loss reduction for a given step is estimated using +!> gradient and hessian +!> \param grad ... +!> \param curvature ... +!> \param trust_step ... +!> \param step_size ... +!> \param border ... +!> \par History +!> 2019.12 created [Rustam Z Khaliullin] +!> \author Rustam Z Khaliullin +! ************************************************************************************************** + SUBROUTINE trust_step_subproblem(grad, curvature, trust_step, step_size, border) + + REAL(KIND=dp), INTENT(IN) :: grad, curvature, trust_step + REAL(KIND=dp), INTENT(INOUT) :: step_size + LOGICAL, INTENT(INOUT) :: border + + REAL(KIND=dp) :: grad_sign + + IF (curvature .LT. 0.0_dp ) THEN + grad_sign = SIGN(1.0_dp, grad) ! sign of grad + step_size=-grad_sign*trust_step + border = .TRUE. + ELSE + step_size = -grad/curvature + IF (step_size .GT. trust_step) THEN + step_size=trust_step + border = .TRUE. + ELSE + border = .FALSE. + ENDIF + ENDIF + + END SUBROUTINE trust_step_subproblem + ! ************************************************************************************************** !> \brief Loss reduction for a given step is estimated using !> gradient and hessian diff --git a/src/nlmo_optimizer.F b/src/nlmo_optimizer.F index d79ef379f12..19653365082 100644 --- a/src/nlmo_optimizer.F +++ b/src/nlmo_optimizer.F @@ -21,26 +21,40 @@ MODULE nlmo_optimizer cp_logger_get_default_unit_nr,& cp_logger_type USE dbcsr_api, ONLY: & - dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_filter, dbcsr_multiply, dbcsr_norm, & - dbcsr_norm_maxabsnorm, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, & - dbcsr_type_no_symmetry - USE input_constants, ONLY: almo_scf_pcg,& - almo_scf_trustr,& - trustr_cauchy,& - trustr_dogleg,& - xalmo_prec_dbfgs,& - xalmo_prec_lbfgs - USE iterate_matrix, ONLY: determinant + dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, & + dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_finalize, & + dbcsr_frobenius_norm, dbcsr_func_dtanh, dbcsr_func_inverse, dbcsr_func_tanh, & + dbcsr_function_of_elements, dbcsr_get_block_p, dbcsr_get_diag, dbcsr_get_info, & + dbcsr_hadamard_product, 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_norm, dbcsr_norm_maxabsnorm, & + dbcsr_p_type, dbcsr_print_block_sum, dbcsr_release, dbcsr_reserve_block2d, dbcsr_scale, & + dbcsr_set, dbcsr_set_diag, dbcsr_triu, dbcsr_type, dbcsr_type_no_symmetry, & + dbcsr_work_create, dbcsr_print + USE input_constants, ONLY: & + 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, cg_bfgs_hybrid, & + 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_identity, & + almo_scf_pcg, almo_scf_skip, almo_scf_trustr, xalmo_prec_lbfgs, xalmo_prec_dbfgs + USE iterate_matrix, ONLY: determinant,& + invert_Hotelling,& + matrix_sqrt_Newton_Schulz USE kinds, ONLY: dp + USE nlmo_types, ONLY: nlmo_env_type USE machine, ONLY: m_walltime USE nlmo_methods, ONLY: & - nlmo_env_allocate_all, nlmo_env_compute_BK_ij, nlmo_env_compute_sigma0_ij, & - nlmo_env_deallocate_all, nlmo_env_hessian_apply, nlmo_env_init, nlmo_env_invert_hessian, & - nlmo_env_loss_curvature_grad_to_step, nlmo_env_loss_function, nlmo_env_loss_gradient, & - nlmo_env_loss_hessian, nlmo_env_mainvar_initial_guess, nlmo_env_mainvar_to_aux, & - nlmo_env_predicted_reduction, nlmo_env_release, nlmo_env_set_flags - USE nlmo_types, ONLY: nlmo_env_type - USE qs_environment_types, ONLY: qs_environment_type + nlmo_env_init, nlmo_env_allocate_all, nlmo_env_compute_BK_ij, nlmo_env_compute_sigma0_ij, & + nlmo_env_deallocate_all, nlmo_env_loss_curvature_grad_to_step, nlmo_env_loss_function, & + nlmo_env_loss_gradient, nlmo_env_mainvar_initial_guess, nlmo_env_mainvar_to_aux, & + nlmo_env_predicted_reduction, nlmo_env_hessian_apply, nlmo_env_minus_hessian_inv_apply, & + nlmo_env_set_flags, nlmo_env_invert_hessian, nlmo_env_loss_hessian, nlmo_env_release, & + trust_step_subproblem + USE qs_energy_types, ONLY: qs_energy_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_loc_utils, ONLY: compute_berry_operator + USE qs_localization_methods, ONLY: initialize_weights #include "./base/base_uses.f90" IMPLICIT NONE @@ -75,11 +89,13 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) routineP = moduleN//':'//routineN CHARACTER(LEN=59) :: print_string - INTEGER :: unit_nr + INTEGER :: unit_nr, handle, iteration REAL(KIND=dp) :: det_diff, penalty_alpha, prev_determinant TYPE(cp_logger_type), POINTER :: logger TYPE(nlmo_env_type) :: nlmo_env + CALL timeset(routineN, handle) + ! get a useful output_unit logger => cp_get_default_logger() IF (logger%para_env%mepos == logger%para_env%source) THEN @@ -88,6 +104,9 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) unit_nr = -1 ENDIF + iteration = 0 + det_diff = 0.0_dp + !create a print function IF (unit_nr > 0) THEN WRITE (unit_nr, *) @@ -136,8 +155,26 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) ! loop over the strength of the orthogonalization penalty prev_determinant = 10.0_dp penalty_alpha = almo_scf_env%opt_nlmo_penalty%penalty_strength + + SELECT CASE (almo_scf_env%construct_nlmos) + CASE (almo_scf_pcg) + IF (almo_scf_env%opt_nlmo_pcg%eps_error_early .NE. -1.0_dp) THEN + nlmo_env%eps_error = almo_scf_env%opt_nlmo_pcg%eps_error_early + ELSE + nlmo_env%eps_error = almo_scf_env%opt_nlmo_pcg%eps_error + ENDIF + CASE (almo_scf_trustr) + IF (almo_scf_env%opt_nlmo_trustr%eps_error_early.NE. -1.0_dp) THEN + nlmo_env%eps_error = almo_scf_env%opt_nlmo_trustr%eps_error_early + ELSE + nlmo_env%eps_error = almo_scf_env%opt_nlmo_trustr%eps_error + ENDIF + END SELECT + DO ! WHILE (almo_scf_env%overlap_determinant .GT. almo_scf_env%opt_nlmo_penalty%final_determinant) + iteration = iteration + 1 + IF (unit_nr > 0) THEN WRITE (unit_nr, '()') WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) @@ -152,6 +189,7 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) CASE (almo_scf_pcg) CALL construct_nlmos_pcg(qs_env=qs_env, & nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & optimizer=almo_scf_env%opt_nlmo_pcg, & matrix_s=almo_scf_env%matrix_s(1), & matrix_mo_inout=almo_scf_env%matrix_t, & @@ -159,6 +197,7 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) CASE (almo_scf_trustr) CALL construct_nlmos_trustr(qs_env=qs_env, & nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & optimizer=almo_scf_env%opt_nlmo_trustr, & matrix_s=almo_scf_env%matrix_s(1), & matrix_mo_inout=almo_scf_env%matrix_t, & @@ -169,6 +208,7 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) CASE (almo_scf_pcg) CALL construct_nlmos_pcg(qs_env=qs_env, & nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & optimizer=almo_scf_env%opt_nlmo_pcg, & matrix_s=almo_scf_env%matrix_s(1), & matrix_mo_inout=almo_scf_env%matrix_v, & @@ -176,6 +216,7 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) CASE (almo_scf_trustr) CALL construct_nlmos_trustr(qs_env=qs_env, & nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & optimizer=almo_scf_env%opt_nlmo_trustr, & matrix_s=almo_scf_env%matrix_s(1), & matrix_mo_inout=almo_scf_env%matrix_v, & @@ -187,9 +228,73 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) almo_scf_env%overlap_determinant = MAXVAL(nlmo_env%overlap_determinant) det_diff = prev_determinant - almo_scf_env%overlap_determinant - IF (det_diff < almo_scf_env%opt_nlmo_penalty%determinant_tolerance .OR. & - almo_scf_env%overlap_determinant .LE. almo_scf_env%opt_nlmo_penalty%final_determinant) THEN + IF ((det_diff .GE. 0.0_dp) .AND. & + (det_diff < almo_scf_env%opt_nlmo_penalty%determinant_tolerance) .OR. & + (almo_scf_env%overlap_determinant .LE. almo_scf_env%opt_nlmo_penalty%final_determinant)) THEN + IF ((almo_scf_env%opt_nlmo_pcg%eps_error_early .NE. -1.0_dp .AND. & + almo_scf_env%opt_nlmo_pcg%eps_error .NE. almo_scf_env%opt_nlmo_pcg%eps_error_early) .OR. & + (almo_scf_env%opt_nlmo_trustr%eps_error_early .NE. -1.0_dp .AND. & + almo_scf_env%opt_nlmo_trustr%eps_error .NE. almo_scf_env%opt_nlmo_trustr%eps_error_early)) THEN + + SELECT CASE (almo_scf_env%construct_nlmos) + CASE (almo_scf_pcg) + nlmo_env%eps_error = almo_scf_env%opt_nlmo_pcg%eps_error + CASE (almo_scf_trustr) + nlmo_env%eps_error = almo_scf_env%opt_nlmo_trustr%eps_error + END SELECT + + IF (unit_nr > 0) THEN + WRITE (unit_nr, '()') + WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) + print_string = "Penalty alpha (dimensionless):" + WRITE (unit_nr, '(T2,A59,F20.10)') print_string, penalty_alpha + WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) + WRITE (unit_nr, '()') + ENDIF + + IF (.NOT. virtuals) THEN + SELECT CASE (almo_scf_env%construct_nlmos) + CASE (almo_scf_pcg) + CALL construct_nlmos_pcg(qs_env=qs_env, & + nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & + optimizer=almo_scf_env%opt_nlmo_pcg, & + matrix_s=almo_scf_env%matrix_s(1), & + matrix_mo_inout=almo_scf_env%matrix_t, & + matrix_templateOO=almo_scf_env%matrix_sigma_inv) + CASE (almo_scf_trustr) + CALL construct_nlmos_trustr(qs_env=qs_env, & + nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & + optimizer=almo_scf_env%opt_nlmo_trustr, & + matrix_s=almo_scf_env%matrix_s(1), & + matrix_mo_inout=almo_scf_env%matrix_t, & + matrix_templateOO=almo_scf_env%matrix_sigma_inv) + END SELECT + ELSE + SELECT CASE (almo_scf_env%construct_nlmos) + CASE (almo_scf_pcg) + CALL construct_nlmos_pcg(qs_env=qs_env, & + nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & + optimizer=almo_scf_env%opt_nlmo_pcg, & + matrix_s=almo_scf_env%matrix_s(1), & + matrix_mo_inout=almo_scf_env%matrix_v, & + matrix_templateOO=almo_scf_env%matrix_sigma_vv) + CASE (almo_scf_trustr) + CALL construct_nlmos_trustr(qs_env=qs_env, & + nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & + optimizer=almo_scf_env%opt_nlmo_trustr, & + matrix_s=almo_scf_env%matrix_s(1), & + matrix_mo_inout=almo_scf_env%matrix_v, & + matrix_templateOO=almo_scf_env%matrix_sigma_vv) + END SELECT + ENDIF + ENDIF + EXIT + ENDIF prev_determinant = almo_scf_env%overlap_determinant @@ -198,6 +303,17 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) penalty_alpha = penalty_alpha/ & ABS(almo_scf_env%opt_nlmo_penalty%penalty_strength_dec_factor) + SELECT CASE (almo_scf_env%construct_nlmos) + CASE (almo_scf_pcg) + IF (almo_scf_env%opt_nlmo_pcg%eps_error_early .NE. -1.0_dp) THEN + nlmo_env%eps_error = nlmo_env%eps_error / 2.0_dp + ENDIF + CASE (almo_scf_trustr) + IF (almo_scf_env%opt_nlmo_trustr%eps_error_early.NE. -1.0_dp) THEN + nlmo_env%eps_error = nlmo_env%eps_error / 2.0_dp + ENDIF + END SELECT + ENDDO IF (unit_nr > 0) THEN @@ -225,12 +341,15 @@ SUBROUTINE nlmo_optimization_entry(qs_env, almo_scf_env, virtuals) CALL nlmo_env_release(nlmo_env) + CALL timestop(handle) + END SUBROUTINE nlmo_optimization_entry ! ************************************************************************************************** !> \brief Optimization of NLMOs using PCG minimizers !> \param qs_env ... !> \param nlmo_env ... +!> \param almo_scf_env ... !> \param optimizer controls the optimization algorithm !> \param matrix_s - AO overlap (NAOs x NAOs) !> \param matrix_mo_inout - initial and final MOs (NAOs x NMOs) @@ -239,11 +358,12 @@ END SUBROUTINE nlmo_optimization_entry !> 2018.10 created [Rustam Z Khaliullin] !> \author Rustam Z Khaliullin ! ************************************************************************************************** - SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & - matrix_s, matrix_mo_inout, & - matrix_templateOO) + SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, almo_scf_env, optimizer, & + matrix_s, matrix_mo_inout, & + matrix_templateOO) TYPE(qs_environment_type), POINTER :: qs_env TYPE(nlmo_env_type), INTENT(INOUT) :: nlmo_env + TYPE(almo_scf_env_type), INTENT(IN) :: almo_scf_env TYPE(optimizer_options_type), INTENT(INOUT) :: optimizer TYPE(dbcsr_type), INTENT(IN) :: matrix_s TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), & @@ -258,20 +378,22 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & INTEGER :: cg_iteration, handle, ispin, iteration, line_search_iteration, & linear_search_type, max_iter, nspins, outer_iteration, outer_max_iter, unit_nr LOGICAL :: converged, just_started, line_search, & - outer_prepare_to_exit, & - prepare_to_exit, reset_conjugator + outer_prepare_to_exit, identity_on, & + prepare_to_exit, reset_conjugator, & + border, rejected REAL(KIND=dp) :: appr_sec_der, beta, denom, denom2, e0, e1, g0, g0sign, g1, g1sign, & - grad_norm, line_search_error, localization_obj_function, next_step_size_guess, & + grad_norm, line_search_error, localization_obj_function, next_step_size_guess, trust_max, & obj_function_ispin, objf_diff, objf_new, objf_old, penalty_func_new, step_size, t1, t2, & - tempreal + tempreal, l_index, trust_step, B0, B1, rho_step, eta_step REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: grad_norm_spin REAL(KIND=dp), DIMENSION(2) :: loss_components, overlap_determinant TYPE(cp_logger_type), POINTER :: logger - TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: prev_grad, prev_m_theta, & + TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: prev_grad, prev_m_theta, Hstep, & prev_minus_prec_grad, prev_step, step CALL timeset(routineN, handle) + ! get a useful output_unit logger => cp_get_default_logger() IF (logger%para_env%mepos == logger%para_env%source) THEN @@ -280,6 +402,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & unit_nr = -1 ENDIF + nspins = SIZE(matrix_mo_inout) CALL nlmo_env_set_flags(nlmo_env=nlmo_env, & optimizer=optimizer) @@ -288,9 +411,11 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & templateNO=matrix_mo_inout) CALL nlmo_env_mainvar_initial_guess(nlmo_env=nlmo_env, & m_mo_initial=matrix_mo_inout) + ! compute the overlap of the initial orbitals CALL nlmo_env_compute_sigma0_ij(nlmo_env=nlmo_env, & matrix_s=matrix_s) + ! use initial orbitals to compute the BK matrix CALL nlmo_env_compute_BK_ij(nlmo_env=nlmo_env, & qs_env=qs_env, & @@ -303,6 +428,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ALLOCATE (prev_step(nspins)) ALLOCATE (step(nspins)) ALLOCATE (prev_minus_prec_grad(nspins)) + ALLOCATE (Hstep(nspins)) DO ispin = 1, nspins @@ -322,6 +448,9 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & CALL dbcsr_create(prev_minus_prec_grad(ispin), & template=nlmo_env%m_sigma0(ispin), & matrix_type=dbcsr_type_no_symmetry) + CALL dbcsr_create(Hstep(ispin), & + template=nlmo_env%m_sigma0(ispin), & + matrix_type=dbcsr_type_no_symmetry) CALL dbcsr_set(step(ispin), 0.0_dp) CALL dbcsr_set(prev_step(ispin), 0.0_dp) @@ -334,8 +463,10 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & 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 + linear_search_type = 3 ! safe restart, no quadratic assumption, takes more steps localization_obj_function = 0.0_dp + identity_on = .TRUE. + l_index = 0.0_dp DO @@ -373,9 +504,9 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ! 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), nlmo_env%grad(ispin)) - ENDIF + IF (line_search_iteration .EQ. 0 .AND. iteration .NE. 0) THEN + CALL dbcsr_copy(prev_grad(ispin), nlmo_env%grad(ispin)) + ENDIF ENDDO ! ispin CALL nlmo_env_loss_gradient(nlmo_env) @@ -387,7 +518,9 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ENDDO ! ispin grad_norm = MAXVAL(grad_norm_spin) - converged = (grad_norm .LE. optimizer%eps_error) + converged = (grad_norm .LE. nlmo_env%eps_error) + + !IF (converged .AND. (.NOT. line_search) .OR. (iteration .GE. max_iter)) THEN IF (converged .OR. (iteration .GE. max_iter)) THEN prepare_to_exit = .TRUE. ENDIF @@ -459,22 +592,36 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & CALL dbcsr_copy(prev_step(ispin), step(ispin)) ENDDO ! ispin + IF (grad_norm .LT. almo_scf_env%opt_nlmo_newton_pcg%start_grad .AND. & + iteration .GT. 100) THEN + identity_on = .FALSE. + ENDIF + + IF (grad_norm .GT. almo_scf_env%opt_nlmo_newton_pcg%start_grad .OR. & + iteration .LT. 100) THEN + identity_on = .TRUE. + ENDIF + ! compute the new step CALL nlmo_env_loss_curvature_grad_to_step(nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & grad=nlmo_env%grad, & step=step, & prev_grad=prev_grad, & prev_m_theta=prev_m_theta, & - iteration=iteration) + iteration=iteration, & + identity_on=identity_on) ! check whether we need to reset conjugate directions IF (iteration .EQ. 0) THEN reset_conjugator = .TRUE. ENDIF ! compute the conjugation coefficient - beta + ! start_grad default value is 0, such that it won't effect IF (.NOT. reset_conjugator) THEN CALL compute_cg_beta( & beta=beta, & + l_index=l_index, & reset_conjugator=reset_conjugator, & conjugator=optimizer%conjugator, & grad=nlmo_env%grad, & @@ -488,6 +635,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & IF (reset_conjugator) THEN beta = 0.0_dp + l_index = 0.0_dp IF (unit_nr > 0 .AND. (.NOT. just_started)) THEN WRITE (unit_nr, '(T2,A35)') "Re-setting conjugator to zero" ENDIF @@ -501,12 +649,18 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & 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) + IF (optimizer%conjugator .EQ. cg_bfgs_hybrid) THEN + CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta) + CALL dbcsr_add(step(ispin), nlmo_env%grad(ispin), 1.0_dp, -1.0_dp*l_index) + ELSE IF (.NOT. reset_conjugator) THEN + CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta) + ENDIF ENDDO ! ispin ENDIF ! update the step direction + ! estimate the step size IF (.NOT. line_search) THEN ! we just changed the direction and @@ -514,12 +668,28 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ! 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(nlmo_env%grad(ispin), step(ispin), tempreal) - g0 = g0 + tempreal - ENDDO ! ispin + B0 = 0.0_dp + IF (linear_search_type .EQ. 4) THEN ! this is trust_step + !CALL nlmo_env_hessian_apply(nlmo_env=nlmo_env, m_in=step, m_out=Hstep) + DO ispin = 1, nspins + CALL dbcsr_dot(nlmo_env%grad(ispin), step(ispin), tempreal) + g0 = g0 + tempreal + !CALL dbcsr_dot(step(ispin), Hstep(ispin), tempreal) + !B0 = B0 + tempreal + ENDDO ! ispin + IF (unit_nr > 0) THEN + !WRITE (unit_nr, '(T2,A59)') "Compute module gradient and curvature" + !WRITE (unit_nr, '(T2,A19,F19.5)') "Curvature:", B0 + WRITE (unit_nr, '(T2,A59)') "The first step use the initial LS guess" + ENDIF + ELSE + DO ispin = 1, nspins + CALL dbcsr_dot(nlmo_env%grad(ispin), step(ispin), tempreal) + g0 = g0 + tempreal + ENDDO ! ispin + ENDIF g0sign = SIGN(1.0_dp, g0) ! sign of g0 - IF (linear_search_type .EQ. 1) THEN ! this is quadratic LS + IF (linear_search_type .EQ. 1 .OR. linear_search_type .EQ. 3) THEN ! this is quadratic LS IF (iteration .EQ. 0) THEN step_size = optimizer%lin_search_step_size_guess ELSE @@ -535,29 +705,132 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ! 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 + ELSE IF (linear_search_type .EQ. 4) THEN ! this is trust_step + trust_max = 2.0_dp + trust_step = MIN(ABS(optimizer%lin_search_step_size_guess),trust_max) + eta_step = 0.25_dp + rejected = .FALSE. + !CALL trust_step_subproblem(grad=g0, & + ! curvature=B0, & + ! trust_step=trust_step, & + ! step_size=step_size, & + ! border=border) + !step_size = optimizer%lin_search_step_size_guess + IF (iteration .EQ. 0) THEN + step_size = 0.000000001_dp + ELSE + IF (next_step_size_guess .LE. 0.0_dp) THEN + step_size = 0.000000001_dp + ELSE + ! take the last value + step_size = optimizer%lin_search_step_size_guess + !step_size = next_step_size_guess*1.05_dp + ENDIF + ENDIF + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A19,F19.8)') "Proposed step size:", step_size + ENDIF + ELSE IF (linear_search_type .EQ. 5) THEN + 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 + !IF (unit_nr > 0) THEN + ! iter_type = TRIM("Line position") + ! WRITE (unit_nr, '(T10,6A19)') iter_type, "Line grad", "Next line step", & + ! "Obj Function", "Lol Function", "Pena Function" + ! iter_type = TRIM("LS:") + ! WRITE (unit_nr, '(T2,A8,6F20.8)') iter_type, 0.0_dp, g0, step_size, & + ! e0, localization_obj_function, penalty_func_new + !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(nlmo_env%grad(ispin), step(ispin), tempreal) - g1 = g1 + tempreal - ENDDO ! ispin + B1 = 0.0_dp + IF (linear_search_type .EQ. 4) THEN ! this is trust_step + !CALL nlmo_env_hessian_apply(nlmo_env=nlmo_env, m_in=step, m_out=Hstep) + DO ispin = 1, nspins + CALL dbcsr_dot(nlmo_env%grad(ispin), step(ispin), tempreal) + g1 = g1 + tempreal + !CALL dbcsr_dot(step(ispin), Hstep(ispin), tempreal) + !B1 = B1 + tempreal + ENDDO ! ispin + B1 = (g1 - g0)/step_size + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A60)') "Compute module gradient and curvature" + WRITE (unit_nr, '(T2,A19,F19.5)') "Curvature:", B1 + ENDIF + ELSE + DO ispin = 1, nspins + CALL dbcsr_dot(nlmo_env%grad(ispin), step(ispin), tempreal) + g1 = g1 + tempreal + ENDDO ! ispin + ENDIF g1sign = SIGN(1.0_dp, g1) ! sign of g1 - IF (linear_search_type .EQ. 1) THEN + IF (linear_search_type .EQ. 4) THEN ! trust_step algo + + IF (rejected) THEN + CALL trust_step_subproblem(grad=g1, & + curvature=B1, & + trust_step=trust_step, & + step_size=step_size, & + border=border) + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A30,F19.5)') "Rejected and re-evaluated step size:", step_size + ENDIF + ELSE + rho_step = ( e1 - e0 ) / ( g0*step_size + 0.5_dp*step_size**2*B0 ) + + IF ( rho_step .LT. 0.25_dp ) THEN + trust_step = 0.25_dp*trust_step + ELSE + IF ( rho_step .GT. 0.75_dp .AND. border) THEN + trust_step = MIN(2.0_dp*trust_step,trust_max) + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A30)') "Step size to border!!!" + ENDIF + ENDIF + ENDIF + + IF ( rho_step .GT. eta_step ) THEN + CALL trust_step_subproblem(grad=g1, & + curvature=B1, & + trust_step=trust_step, & + step_size=step_size, & + border=border) + rejected = .FALSE. + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A60)') "Accepted" + ENDIF + ELSE + step_size = -step_size + border = .FALSE. + rejected = .TRUE. + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2, A30)') "Rejected and Redo" + ENDIF + ENDIF + ENDIF + + ELSE 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 (g1sign .NE. g0sign) THEN + step_size = -step_size/2.0; + ELSE + !appr_sec_der = (g1 - g0)/step_size + IF (appr_sec_der .GT. 0.0_dp) THEN + step_size = -g1/appr_sec_der + ELSE + step_size = step_size*1.5; + ENDIF + ENDIF !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 + ! WRITE (unit_nr, '(A2,T5,F30.5, 2F19.5)') & + ! "DT", appr_sec_der, step_size, -g1/appr_sec_der + ! !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 @@ -566,12 +839,26 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ELSE step_size = step_size*1.5; ENDIF + ELSE IF (linear_search_type .EQ. 3) THEN + 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. 5) THEN + step_size = optimizer%lin_search_step_size_guess 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 + !IF (unit_nr > 0) THEN + ! iter_type = TRIM("Line position") + ! WRITE (unit_nr, '(T10,6A19)') iter_type, "Line grad", "Next line step", & + ! "Obj Function", "Lol Function", "Pena Function" + ! iter_type = TRIM("LS:") + ! WRITE (unit_nr, '(T2,A8,6F20.8)') iter_type, next_step_size_guess, g1, step_size, & + ! e1, localization_obj_function, penalty_func_new + !ENDIF e0 = e1 g0 = g1 g0sign = g1sign @@ -580,6 +867,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & ! update theta DO ispin = 1, nspins + CALL dbcsr_print(nlmo_env%m_theta(ispin)) 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), nlmo_env%m_theta(ispin)) @@ -598,7 +886,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & 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)') & + WRITE (unit_nr, '(T2,A13,I6,F23.10,E14.5,ES14.5,F10.5)') & iter_type, iteration, & objf_new, objf_diff, grad_norm, & t2 - t1 @@ -630,6 +918,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & CALL dbcsr_release(prev_step(ispin)) CALL dbcsr_release(step(ispin)) CALL dbcsr_release(prev_minus_prec_grad(ispin)) + CALL dbcsr_release(Hstep(ispin)) ENDDO ! ispin DEALLOCATE (grad_norm_spin) @@ -638,6 +927,7 @@ SUBROUTINE construct_nlmos_pcg(qs_env, nlmo_env, optimizer, & DEALLOCATE (prev_step) DEALLOCATE (step) DEALLOCATE (prev_minus_prec_grad) + DEALLOCATE (Hstep) DO ispin = 1, nspins CALL dbcsr_copy(matrix_mo_inout(ispin), nlmo_env%m_mo(ispin)) @@ -657,6 +947,7 @@ END SUBROUTINE construct_nlmos_pcg !> \brief Optimization of ALMOs using trust region minimizers !> \param qs_env ... !> \param nlmo_env ... +!> \param almo_scf_env ... !> \param optimizer controls the optimization algorithm !> \param matrix_s ... !> \param matrix_mo_inout ... @@ -665,12 +956,13 @@ END SUBROUTINE construct_nlmos_pcg !> 2020.01 created [Rustam Z Khaliullin] !> \author Rustam Z Khaliullin ! ************************************************************************************************** - SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & + SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, almo_scf_env, optimizer, & matrix_s, matrix_mo_inout, & matrix_templateOO) TYPE(qs_environment_type), POINTER :: qs_env TYPE(nlmo_env_type), INTENT(INOUT) :: nlmo_env + TYPE(almo_scf_env_type), INTENT(IN) :: almo_scf_env TYPE(optimizer_options_type), INTENT(INOUT) :: optimizer TYPE(dbcsr_type), INTENT(IN) :: matrix_s TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), & @@ -686,12 +978,12 @@ SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & outer_iteration, unit_nr LOGICAL :: border_reached, inner_loop_success, & reset_conjugator, same_position, & - scf_converged + scf_converged, identity_on REAL(kind=dp) :: beta, energy_start, energy_trial, eta, expected_reduction, & fake_step_size_to_report, grad_norm_ratio, grad_norm_ref, loss_change_to_report, & - loss_start, loss_trial, model_grad_norm, penalty_start, penalty_trial, radius_current, & - radius_max, real_temp, rho, spin_factor, step_norm, step_size, t1, t1outer, t2, t2outer, & - y_scalar + loss_start, loss_trial, model_grad_norm, penalty_start, penalty_trial, & + radius_current, radius_max, real_temp, rho, spin_factor, step_norm, step_size, t1, & + t1outer, t2, t2outer, y_scalar REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: grad_norm_spin REAL(KIND=dp), DIMENSION(2) :: loss_components, overlap_determinant TYPE(cp_logger_type), POINTER :: logger @@ -712,6 +1004,7 @@ SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & unit_nr = -1 ENDIF + nspins = SIZE(matrix_mo_inout) CALL nlmo_env_set_flags(nlmo_env=nlmo_env, & optimizer=optimizer) @@ -798,6 +1091,7 @@ SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & penalty_trial = 0.0_dp loss_start = 0.0_dp ! sum of the energy and penalty loss_trial = 0.0_dp + identity_on = .FALSE. same_position = .FALSE. @@ -849,9 +1143,16 @@ SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & radius=radius_current, & new=.NOT. same_position, & time=t2outer - t1outer) + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A19,F23.10)') & + "Localization:", energy_start + WRITE (unit_nr, '(T2,A19,F23.10)') & + "Orthogonalization:", penalty_start + ENDIF + t1outer = m_walltime() - IF (grad_norm_ref .LE. optimizer%eps_error) THEN + IF (grad_norm_ref .LE. nlmo_env%eps_error) THEN scf_converged = .TRUE. border_reached = .FALSE. expected_reduction = 0.0_dp @@ -1026,7 +1327,10 @@ SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & ENDIF - IF (optimizer%trustr_algorithm .EQ. trustr_cauchy) THEN + !IF (optimizer%trustr_algorithm .EQ. trustr_cauchy) THEN + IF (optimizer%trustr_algorithm .EQ. trustr_cauchy .OR. & + (optimizer%trustr_algorithm .EQ. trustr_dogleg .AND. & + grad_norm_ref .GT. almo_scf_env%opt_nlmo_newton_pcg%start_grad)) THEN ! trustr_steihaug, trustr_cauchy, trustr_dogleg border_reached = .FALSE. @@ -1052,19 +1356,33 @@ SUBROUTINE construct_nlmos_trustr(qs_env, nlmo_env, optimizer, & ! invert B IF (unit_nr > 0 .AND. debug_mode) WRITE (unit_nr, *) "...(Pseudo-)invert model Hessian" - CALL nlmo_env_invert_hessian(nlmo_env) - - ! get pB = Binv.m_model_r = -Binv.grad - ! RZK-critical: implement an nlmo_env method - DO ispin = 1, nspins - - CALL dbcsr_multiply("N", "N", 1.0_dp, & - nlmo_env%m_model_hessian_inv(ispin), & - m_model_r(ispin), & - 0.0_dp, m_model_Bd(ispin), & - filter_eps=nlmo_env%eps_filter) - - ENDDO ! ispin + !CALL nlmo_env_invert_hessian(nlmo_env) + + !! get pB = Binv.m_model_r = -Binv.grad + !! RZK-critical: implement an nlmo_env method + !DO ispin = 1, nspins + + ! IF (iteration .GT. 0) THEN + ! CALL dbcsr_scale(nlmo_env%m_model_hessian_inv(ispin), -1.0_dp) + ! END IF + + ! CALL dbcsr_multiply("N", "N", 1.0_dp, & + ! nlmo_env%m_model_hessian_inv(ispin), & + ! m_model_r(ispin), & + ! 0.0_dp, m_model_Bd(ispin), & + ! filter_eps=nlmo_env%eps_filter) + + !ENDDO ! ispin + + ! Compute pB through B * pB = -grad by conjugate gradient methods. + CALL nlmo_env_loss_curvature_grad_to_step(nlmo_env=nlmo_env, & + almo_scf_env=almo_scf_env, & + grad=nlmo_env%grad, & + step=m_model_Bd, & + prev_grad=m_model_r_prev, & ! not important here + prev_m_theta=prev_step, & ! not important here + iteration=iteration, & + identity_on=identity_on) ! Compute norm of pB CALL contravariant_matrix_norm( & diff --git a/src/nlmo_types.F b/src/nlmo_types.F index f05346fdd96..80c9bbccab9 100644 --- a/src/nlmo_types.F +++ b/src/nlmo_types.F @@ -29,6 +29,7 @@ MODULE nlmo_types INTEGER :: loc_operator, & natoms REAL(KIND=dp) :: eps_filter + REAL(KIND=dp) :: eps_error ! AO matrices TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set diff --git a/src/qs_localization_methods.F b/src/qs_localization_methods.F index 896dc5dd59f..5a1a190fcc1 100644 --- a/src/qs_localization_methods.F +++ b/src/qs_localization_methods.F @@ -1291,6 +1291,12 @@ SUBROUTINE direct_mini(weights, zij, vectors, max_iter, eps_localization, iterat ENDIF ! lsr .eq. 0 ENDIF ! first step END SELECT + + !IF (output_unit > 0) THEN + ! WRITE (output_unit, '(T5,I10,T18,I10,T31,2F20.6,F10.3)') line_searches, Iterations, Omega, tol, ds_min + ! CALL m_flush(output_unit) + !ENDIF + ! now go to the suggested point ds_min = pos(line_search_count + 1) ds = pos(line_search_count + 1) - pos(line_search_count) diff --git a/tests/QS/regtest-nlmo/Si-nlmos.inp b/tests/QS/regtest-nlmo/Si-nlmos.inp index fd964e207b0..3126ff20773 100644 --- a/tests/QS/regtest-nlmo/Si-nlmos.inp +++ b/tests/QS/regtest-nlmo/Si-nlmos.inp @@ -31,6 +31,7 @@ XALMO_TRIAL_WF SIMPLE CONSTRUCT_NLMOS PCG + OCCUPIED_NLMOS TRUE VIRTUAL_NLMOS FALSE NLMO_COMPACT_FILTER_START 1.0E-2 NLMO_OPERATOR PIPEK @@ -57,7 +58,7 @@ &NLMO_PENALTY PENALTY_STRENGTH 0.01 - DETERMINANT_TOLERANCE 1.0E-8 + DETERMINANT_TOLERANCE 1.0E-4 PENALTY_STRENGTH_DECREASE_FACTOR 2.0 FINAL_DETERMINANT 0.6 &END NLMO_PENALTY diff --git a/tests/QS/regtest-nlmo/TEST_FILES b/tests/QS/regtest-nlmo/TEST_FILES index 0743c19dd5d..d12117e8e94 100644 --- a/tests/QS/regtest-nlmo/TEST_FILES +++ b/tests/QS/regtest-nlmo/TEST_FILES @@ -1,4 +1,9 @@ -pipek_C6H6.inp 90 1e-04 169.9680388838 +pipek_C6H6.inp 90 1e-04 172.6818129438 +pipek_C6H6_CG_BFGS.inp 90 1e-04 172.6859039706 +pipek_C6H6_trustr.inp 90 1e-04 172.6828872753 Si-nlmos.inp 90 1e-04 118.5296581594 water_berry.inp 90 1e-04 224.3891398828 +water_berry_CG_BFGS.inp 90 1e-04 224.3891388599 +water_berry_trustr.inp 90 1e-04 224.3890270442 +water_berry_trustrDL.inp 90 1e-04 224.3893450367 #EOF diff --git a/tests/QS/regtest-nlmo/pipek_C6H6.inp b/tests/QS/regtest-nlmo/pipek_C6H6.inp index c79194aca8e..6c088e4832e 100644 --- a/tests/QS/regtest-nlmo/pipek_C6H6.inp +++ b/tests/QS/regtest-nlmo/pipek_C6H6.inp @@ -28,6 +28,7 @@ RETURN_ORTHOGONALIZED_MOS FALSE CONSTRUCT_NLMOS PCG + OCCUPIED_NLMOS FALSE VIRTUAL_NLMOS TRUE !FALSE NLMO_OPERATOR PIPEK !NLMO_COMPACT_FILTER_START 1.0E-3 diff --git a/tests/QS/regtest-nlmo/pipek_C6H6_CG_BFGS.inp b/tests/QS/regtest-nlmo/pipek_C6H6_CG_BFGS.inp new file mode 100644 index 00000000000..3ae3be56823 --- /dev/null +++ b/tests/QS/regtest-nlmo/pipek_C6H6_CG_BFGS.inp @@ -0,0 +1,105 @@ +&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 + RETURN_ORTHOGONALIZED_MOS FALSE + + CONSTRUCT_NLMOS PCG + OCCUPIED_NLMOS FALSE + VIRTUAL_NLMOS TRUE !FALSE + NLMO_OPERATOR PIPEK + !NLMO_COMPACT_FILTER_START 1.0E-3 + + &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 200 + EPS_ERROR 0.05 + CONJUGATOR BFGS_HYBRID + PRECONDITIONER LBFGS + LIN_SEARCH_EPS_ERROR 0.1 + LIN_SEARCH_STEP_SIZE_GUESS 0.00001 + MAX_ITER_OUTER_LOOP 5 + &END NLMO_OPTIMIZER_PCG + + &NLMO_PENALTY + PENALTY_STRENGTH 0.05 + DETERMINANT_TOLERANCE 1.0E-8 + PENALTY_STRENGTH_DECREASE_FACTOR 3.0 + FINAL_DETERMINANT 0.8 + &END NLMO_PENALTY + + &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/QS/regtest-nlmo/pipek_C6H6_trustr.inp b/tests/QS/regtest-nlmo/pipek_C6H6_trustr.inp index 4a23d333f55..f99605a80b8 100644 --- a/tests/QS/regtest-nlmo/pipek_C6H6_trustr.inp +++ b/tests/QS/regtest-nlmo/pipek_C6H6_trustr.inp @@ -28,6 +28,7 @@ RETURN_ORTHOGONALIZED_MOS FALSE CONSTRUCT_NLMOS TRUST_REGION + OCCUPIED_NLMOS FALSE VIRTUAL_NLMOS TRUE !FALSE NLMO_OPERATOR PIPEK !NLMO_COMPACT_FILTER_START 1.0E-3 @@ -56,8 +57,8 @@ ALGORITHM CG MAX_ITER 100 MAX_ITER_OUTER_LOOP 1000 - EPS_ERROR 0.00001 - MODEL_GRAD_NORM_RATIO 0.1 + EPS_ERROR 0.05 + MODEL_GRAD_NORM_RATIO 0.7 CONJUGATOR FLETCHER PRECONDITIONER FULL ETA 0.25 @@ -69,7 +70,7 @@ PENALTY_STRENGTH 0.05 DETERMINANT_TOLERANCE 1.0E-8 PENALTY_STRENGTH_DECREASE_FACTOR 3.0 - FINAL_DETERMINANT 0.8 + FINAL_DETERMINANT 0.9 &END NLMO_PENALTY &END ALMO_SCF diff --git a/tests/QS/regtest-nlmo/water_berry.inp b/tests/QS/regtest-nlmo/water_berry.inp index f10703bc73d..c88dfb614e9 100644 --- a/tests/QS/regtest-nlmo/water_berry.inp +++ b/tests/QS/regtest-nlmo/water_berry.inp @@ -31,6 +31,7 @@ XALMO_TRIAL_WF SIMPLE CONSTRUCT_NLMOS PCG + OCCUPIED_NLMOS TRUE VIRTUAL_NLMOS FALSE !TRUE NLMO_COMPACT_FILTER_START 1.0E-2 NLMO_OPERATOR BERRY diff --git a/tests/QS/regtest-nlmo/water_berry_CG_BFGS.inp b/tests/QS/regtest-nlmo/water_berry_CG_BFGS.inp new file mode 100644 index 00000000000..391c371ce5a --- /dev/null +++ b/tests/QS/regtest-nlmo/water_berry_CG_BFGS.inp @@ -0,0 +1,95 @@ +&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 PCG + OCCUPIED_NLMOS TRUE + VIRTUAL_NLMOS FALSE !TRUE + NLMO_COMPACT_FILTER_START 1.0E-2 + NLMO_OPERATOR BERRY + + &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.001 + CONJUGATOR BFGS_HYBRID + PRECONDITIONER LBFGS + LIN_SEARCH_EPS_ERROR 0.5 + LIN_SEARCH_STEP_SIZE_GUESS 0.001 + MAX_ITER_OUTER_LOOP 0 + &END NLMO_OPTIMIZER_PCG + + &NLMO_PENALTY + PENALTY_STRENGTH 0.4 + DETERMINANT_TOLERANCE 1.0E-8 + PENALTY_STRENGTH_DECREASE_FACTOR 2.0 + FINAL_DETERMINANT 0.9 + &END NLMO_PENALTY + + &END ALMO_SCF + + &XC + &XC_FUNCTIONAL BLYP + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 H2O + H 0.000000 -0.757136 0.520545 H2O + H 0.000000 0.757136 0.520545 H2O + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q6 + &END KIND + &TOPOLOGY + CONNECTIVITY GENERATE + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL + diff --git a/tests/QS/regtest-nlmo/water_berry_trustr.inp b/tests/QS/regtest-nlmo/water_berry_trustr.inp index 53906e911d0..d2669c16956 100644 --- a/tests/QS/regtest-nlmo/water_berry_trustr.inp +++ b/tests/QS/regtest-nlmo/water_berry_trustr.inp @@ -31,6 +31,7 @@ XALMO_TRIAL_WF SIMPLE CONSTRUCT_NLMOS TRUST_REGION + OCCUPIED_NLMOS TRUE VIRTUAL_NLMOS FALSE !TRUE NLMO_COMPACT_FILTER_START 1.0E-2 NLMO_OPERATOR BERRY diff --git a/tests/QS/regtest-nlmo/water_berry_trustrDL.inp b/tests/QS/regtest-nlmo/water_berry_trustrDL.inp new file mode 100644 index 00000000000..99e39a3e3db --- /dev/null +++ b/tests/QS/regtest-nlmo/water_berry_trustrDL.inp @@ -0,0 +1,104 @@ +&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 TRUST_REGION + OCCUPIED_NLMOS TRUE + VIRTUAL_NLMOS FALSE !TRUE + NLMO_COMPACT_FILTER_START 1.0E-2 + NLMO_OPERATOR BERRY + + &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_TRUSTR + ALGORITHM Dogleg + MAX_ITER_OUTER_LOOP 1000 + EPS_ERROR 0.001 + CONJUGATOR FLETCHER + PRECONDITIONER FULL + ETA 0.25 + INITIAL_TRUST_RADIUS 0.01 + MAX_TRUST_RADIUS 2.0 + &END NLMO_OPTIMIZER_TRUSTR + + &NLMO_PENALTY + PENALTY_STRENGTH 0.4 + DETERMINANT_TOLERANCE 1.0E-8 + PENALTY_STRENGTH_DECREASE_FACTOR 2.0 + FINAL_DETERMINANT 0.9 + &END NLMO_PENALTY + + &NLMO_NEWTON_PCG + PRECONDITIONER NONE + MAX_ITER_OUTER_LOOP 1 + MAX_ITER 2000 + EPS_ERROR 1.0E-1 + START_GRADIENT 0.01 + &END NLMO_NEWTON_PCG + + &END ALMO_SCF + + &XC + &XC_FUNCTIONAL BLYP + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 H2O + H 0.000000 -0.757136 0.520545 H2O + H 0.000000 0.757136 0.520545 H2O + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q6 + &END KIND + &TOPOLOGY + CONNECTIVITY GENERATE + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL +