Skip to content

Commit

Permalink
Add trust region optimization of block-diagonal ALMOs
Browse files Browse the repository at this point in the history 
Correct previous commit: Add trust region optimization of block-diagonal ALMOs

Detail a dogleg error message
  • Loading branch information
@oschuett
Rustam Z. Khaliullin authored and oschuett committed Jun 20, 2020
1 parent 7cdfc81 commit 06b0715
Show file tree
Hide file tree
Showing 8 changed files with 178 additions and 42 deletions.
59 changes: 31 additions & 28 deletions src/almo_scf.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,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_block_diag_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
Expand Down Expand Up @@ -949,6 +950,9 @@ SUBROUTINE almo_scf_print_job_info(almo_scf_env, unit_nr)
CASE (almo_scf_pcg)
! print out PCG options
CALL print_optimizer_options(almo_scf_env%opt_block_diag_pcg, unit_nr)
CASE (almo_scf_trustr)
! print out TRUST REGION options
CALL print_optimizer_options(almo_scf_env%opt_block_diag_trustr, unit_nr)
END SELECT
ENDIF
Expand Down Expand Up @@ -1288,18 +1292,34 @@ SUBROUTINE almo_scf_main(qs_env, almo_scf_env)
ENDIF

SELECT CASE (almo_scf_env%almo_update_algorithm)
CASE (almo_scf_pcg)
CASE (almo_scf_pcg, almo_scf_trustr, almo_scf_skip)

! ALMO PCG optimizer as a special case of XALMO PCG
CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
almo_scf_env=almo_scf_env, &
optimizer=almo_scf_env%opt_block_diag_pcg, &
quench_t=almo_scf_env%quench_t_blk, &
matrix_t_in=almo_scf_env%matrix_t_blk, &
matrix_t_out=almo_scf_env%matrix_t_blk, &
assume_t0_q0x=.FALSE., &
perturbation_only=.FALSE., &
special_case=xalmo_case_block_diag)
SELECT CASE (almo_scf_env%almo_update_algorithm)
CASE (almo_scf_pcg)

! ALMO PCG optimizer as a special case of XALMO PCG
CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
almo_scf_env=almo_scf_env, &
optimizer=almo_scf_env%opt_block_diag_pcg, &
quench_t=almo_scf_env%quench_t_blk, &
matrix_t_in=almo_scf_env%matrix_t_blk, &
matrix_t_out=almo_scf_env%matrix_t_blk, &
assume_t0_q0x=.FALSE., &
perturbation_only=.FALSE., &
special_case=xalmo_case_block_diag)

CASE (almo_scf_trustr)

CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
almo_scf_env=almo_scf_env, &
optimizer=almo_scf_env%opt_block_diag_trustr, &
quench_t=almo_scf_env%quench_t_blk, &
matrix_t_in=almo_scf_env%matrix_t_blk, &
matrix_t_out=almo_scf_env%matrix_t_blk, &
perturbation_only=.FALSE., &
special_case=xalmo_case_block_diag)

END SELECT

DO ispin = 1, almo_scf_env%nspins
CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
Expand All @@ -1313,30 +1333,13 @@ SUBROUTINE almo_scf_main(qs_env, almo_scf_env)
eps_lanczos=almo_scf_env%eps_lanczos, &
max_iter_lanczos=almo_scf_env%max_iter_lanczos)
ENDDO
!CALL almo_scf_t_blk_to_t_blk_orthonormal(almo_scf_env)

CASE (almo_scf_diag)

! mixing/DIIS optimizer
CALL almo_scf_block_diagonal(qs_env, almo_scf_env, &
almo_scf_env%opt_block_diag_diis)

CASE (almo_scf_skip)

DO ispin = 1, almo_scf_env%nspins
CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
overlap=almo_scf_env%matrix_sigma_blk(ispin), &
metric=almo_scf_env%matrix_s_blk(1), &
retain_locality=.TRUE., &
only_normalize=.FALSE., &
nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
eps_filter=almo_scf_env%eps_filter, &
order_lanczos=almo_scf_env%order_lanczos, &
eps_lanczos=almo_scf_env%eps_lanczos, &
max_iter_lanczos=almo_scf_env%max_iter_lanczos)
ENDDO
!CALL almo_scf_t_blk_to_t_blk_orthonormal(almo_scf_env)

END SELECT

! we might need a copy of the converged KS and sigma_inv
Expand Down
33 changes: 31 additions & 2 deletions src/almo_scf_env_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,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, &
penalty_section, xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, &
almo_opt_pcg_section, almo_opt_trustr_section, almo_scf_section, matrix_iterate_section, &
nlmo_opt_pcg_section, penalty_section, xalmo_opt_newton_pcg_section, xalmo_opt_pcg_section, &
xalmo_opt_trustr_section

CALL timeset(routineN, handle)
Expand All @@ -117,6 +117,8 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env)
"ALMO_OPTIMIZER_DIIS")
almo_opt_pcg_section => section_vals_get_subs_vals(almo_scf_section, &
"ALMO_OPTIMIZER_PCG")
almo_opt_trustr_section => section_vals_get_subs_vals(almo_scf_section, &
"ALMO_OPTIMIZER_TRUSTR")
xalmo_opt_pcg_section => section_vals_get_subs_vals(almo_scf_section, &
"XALMO_OPTIMIZER_PCG")
xalmo_opt_trustr_section => section_vals_get_subs_vals(almo_scf_section, &
Expand Down Expand Up @@ -198,6 +200,33 @@ SUBROUTINE almo_scf_init_read_write_input(input, almo_scf_env)
CALL section_vals_val_get(almo_opt_pcg_section, "PRECONDITIONER", &
i_val=almo_scf_env%opt_block_diag_pcg%preconditioner)

CALL section_vals_val_get(almo_opt_trustr_section, "EPS_ERROR", &
r_val=almo_scf_env%opt_block_diag_trustr%eps_error)
CALL section_vals_val_get(almo_opt_trustr_section, "MAX_ITER", &
i_val=almo_scf_env%opt_block_diag_trustr%max_iter)
CALL section_vals_val_get(almo_opt_trustr_section, "ALGORITHM", &
i_val=almo_scf_env%opt_block_diag_trustr%trustr_algorithm)
CALL section_vals_val_get(almo_opt_trustr_section, "EPS_ERROR_EARLY", &
r_val=almo_scf_env%opt_block_diag_trustr%eps_error_early)
CALL section_vals_val_get(almo_opt_trustr_section, "MAX_ITER_EARLY", &
i_val=almo_scf_env%opt_block_diag_trustr%max_iter_early)
CALL section_vals_val_get(almo_opt_trustr_section, "MAX_ITER_OUTER_LOOP", &
i_val=almo_scf_env%opt_block_diag_trustr%max_iter_outer_loop)
CALL section_vals_val_get(almo_opt_trustr_section, "ETA", &
r_val=almo_scf_env%opt_block_diag_trustr%rho_do_not_update)
almo_scf_env%opt_block_diag_trustr%rho_do_not_update = &
MIN(MAX(almo_scf_env%opt_block_diag_trustr%rho_do_not_update, 0.0_dp), 0.25_dp)
CALL section_vals_val_get(almo_opt_trustr_section, "MODEL_GRAD_NORM_RATIO", &
r_val=almo_scf_env%opt_block_diag_trustr%model_grad_norm_ratio)
CALL section_vals_val_get(almo_opt_trustr_section, "INITIAL_TRUST_RADIUS", &
r_val=almo_scf_env%opt_block_diag_trustr%initial_trust_radius)
CALL section_vals_val_get(almo_opt_trustr_section, "MAX_TRUST_RADIUS", &
r_val=almo_scf_env%opt_block_diag_trustr%max_trust_radius)
CALL section_vals_val_get(almo_opt_trustr_section, "CONJUGATOR", &
i_val=almo_scf_env%opt_block_diag_trustr%conjugator)
CALL section_vals_val_get(almo_opt_trustr_section, "PRECONDITIONER", &
i_val=almo_scf_env%opt_block_diag_trustr%preconditioner)

CALL section_vals_val_get(xalmo_opt_trustr_section, "EPS_ERROR", &
r_val=almo_scf_env%opt_xalmo_trustr%eps_error)
CALL section_vals_val_get(xalmo_opt_trustr_section, "MAX_ITER", &
Expand Down
12 changes: 9 additions & 3 deletions src/almo_scf_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2234,9 +2234,15 @@ END SUBROUTINE apply_domain_operators
!> 2013.01 created [Rustam Z. Khaliullin]
!> \author Rustam Z. Khaliullin
! **************************************************************************************************
SUBROUTINE construct_domain_preconditioner(matrix_main, subm_s_inv, subm_s_inv_half, &
subm_s_half, subm_r_down, matrix_trimmer, dpattern, map, node_of_domain, preconditioner, &
bad_modes_projector_down, use_trimmer, eps_zero_eigenvalues, my_action, skip_inversion)
SUBROUTINE construct_domain_preconditioner(matrix_main, subm_s_inv, &
subm_s_inv_half, subm_s_half, &
subm_r_down, matrix_trimmer, &
dpattern, map, node_of_domain, &
preconditioner, &
bad_modes_projector_down, &
use_trimmer, &
eps_zero_eigenvalues, &
my_action, skip_inversion)
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_main
TYPE(domain_submatrix_type), DIMENSION(:), &
Expand Down
25 changes: 20 additions & 5 deletions src/almo_scf_optimizer.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7001,7 +7001,8 @@ SUBROUTINE compute_preconditioner(domain_prec_out, m_prec_out, m_ks, m_s, &
CHARACTER(len=*), PARAMETER :: routineN = 'compute_preconditioner', &
routineP = moduleN//':'//routineN
INTEGER :: handle, precond_domain_projector
INTEGER :: handle, ndim, precond_domain_projector
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: nn_diagonal
TYPE(dbcsr_type) :: m_tmp_nn_1, m_tmp_no_3
CALL timeset(routineN, handle)
Expand Down Expand Up @@ -7101,14 +7102,26 @@ SUBROUTINE compute_preconditioner(domain_prec_out, m_prec_out, m_ks, m_s, &
IF (special_case .EQ. xalmo_case_block_diag) THEN ! non-overlapping diagonal blocks
IF (skip_inversion) THEN
CPABORT("NYI: impose blk structure on m_prec_out first")
! impose block-diagonal structure
CALL dbcsr_get_info(m_s, nfullrows_total=ndim)
ALLOCATE (nn_diagonal(ndim))
CALL dbcsr_get_diag(m_s, nn_diagonal)
CALL dbcsr_set(m_prec_out, 0.0_dp)
CALL dbcsr_set_diag(m_prec_out, nn_diagonal)
CALL dbcsr_filter(m_prec_out, eps_filter)
DEALLOCATE (nn_diagonal)
CALL dbcsr_copy(m_prec_out, m_tmp_nn_1, keep_sparsity=.TRUE.)
ELSE
CALL pseudo_invert_diagonal_blk( &
matrix_in=m_tmp_nn_1, &
matrix_out=m_prec_out, &
nocc=nocc_of_domain(:) &
)
ENDIF
ELSE IF (special_case .EQ. xalmo_case_fully_deloc) THEN ! the entire system is a block
Expand Down Expand Up @@ -9918,8 +9931,10 @@ SUBROUTINE almo_scf_xalmo_trustr(qs_env, almo_scf_env, optimizer, quench_t, &
eps_filter_in=almo_scf_env%eps_filter &
)
IF (unit_nr > 0 .AND. debug_mode) WRITE (unit_nr, *) "...Step size to border: ", step_size
IF (step_size .GE. 1.0_dp .OR. step_size .LE. 0.0_dp) THEN
CPABORT("Wrong dog-leg step, should not be here")
IF (step_size .GT. 1.0_dp .OR. step_size .LT. 0.0_dp) THEN
IF (unit_nr > 0) &
WRITE (unit_nr, *) "Step size (", step_size, ") must lie inside (0,1)"
CPABORT("Wrong dog leg step. We should never end up here.")
ENDIF
DO ispin = 1, nspins
Expand Down Expand Up @@ -10699,7 +10714,7 @@ SUBROUTINE predicted_reduction(reduction_out, grad_in, step_in, hess_in, &
ENDDO ! ispin
!RZK-critical: do we need to multiply by the spin factor?
!my_reduction = spin_factor*my_reduction
my_reduction = spin_factor*my_reduction
reduction_out = my_reduction
Expand Down
5 changes: 5 additions & 0 deletions src/almo_scf_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -369,6 +369,7 @@ MODULE almo_scf_types
TYPE(optimizer_options_type) :: opt_xalmo_pcg
TYPE(optimizer_options_type) :: opt_xalmo_trustr
TYPE(optimizer_options_type) :: opt_nlmo_pcg
TYPE(optimizer_options_type) :: opt_block_diag_trustr
TYPE(optimizer_options_type) :: opt_xalmo_newton_pcg_solver
TYPE(optimizer_options_type) :: opt_k_pcg

Expand Down Expand Up @@ -504,6 +505,10 @@ SUBROUTINE print_optimizer_options(optimizer, unit_nr)
optimizer%lin_search_step_size_guess
WRITE (unit_nr, '(T4,A,T48,E33.3)') "line search target error:", &
optimizer%lin_search_eps_error
IF (optimizer%neglect_threshold .GT. 0.0_dp) THEN
WRITE (unit_nr, '(T4,A,T48,E33.3)') "low-curvature threshold:", &
optimizer%neglect_threshold
ENDIF

ENDIF

Expand Down
21 changes: 17 additions & 4 deletions src/input_cp2k_almo.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ MODULE input_cp2k_almo
INTEGER, PARAMETER, PRIVATE :: optimizer_xalmo_trustr = 4
INTEGER, PARAMETER, PRIVATE :: optimizer_newton_pcg_solver = 5
INTEGER, PARAMETER, PRIVATE :: optimizer_nlmo_pcg = 6
INTEGER, PARAMETER, PRIVATE :: optimizer_block_diagonal_trustr = 7

PUBLIC :: create_almo_scf_section

Expand Down Expand Up @@ -135,16 +136,17 @@ SUBROUTINE create_almo_scf_section(section)
usage="ALMO_ALGORITHM DIAG", &
default_i_val=almo_scf_diag, &
!enum_c_vals=s2a("DIAG", "DM_SIGN","PCG"),&
!enum_c_vals=s2a("DIAG", "PCG"), &
enum_c_vals=s2a("DIAG", "PCG", "SKIP"), &
enum_c_vals=s2a("DIAG", "PCG", "TRUST_REGION", "SKIP"), &
enum_desc=s2a("DIIS-accelerated diagonalization controlled by ALMO_OPTIMIZER_DIIS. "// &
"Recommended for large systems containing small fragments.", &
!"Update the density matrix using linear scaling routines. "//&
!"Recommended if large fragments are present.",&
"Energy minimization with a PCG algorithm controlled by ALMO_OPTIMIZER_PCG.", &
"Trust-region optimizer is recommended if PCG has difficulty converging. "// &
"It is controlled by ALMO_OPTIMIZER_TRUSTR", &
"Skip optimization of block-diagonal ALMOs."), &
!enum_i_vals=(/almo_scf_diag,almo_scf_dm_sign,almo_scf_pcg/),&
enum_i_vals=(/almo_scf_diag, almo_scf_pcg, almo_scf_skip/))
enum_i_vals=(/almo_scf_diag, almo_scf_pcg, almo_scf_trustr, almo_scf_skip/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

Expand Down Expand Up @@ -529,6 +531,11 @@ SUBROUTINE create_almo_scf_section(section)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

NULLIFY (subsection)
CALL create_optimizer_section(subsection, optimizer_block_diagonal_trustr)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

NULLIFY (subsection)
CALL create_optimizer_section(subsection, optimizer_xalmo_pcg)
CALL section_add_subsection(section, subsection)
Expand Down Expand Up @@ -620,6 +627,12 @@ RECURSIVE SUBROUTINE create_optimizer_section(section, optimizer_id)
"used and controlled by the inner loop", &
n_keywords=10, n_subsections=0, repeats=.FALSE.)
optimizer_type = optimizer_trustr
CASE (optimizer_block_diagonal_trustr)
CALL section_create(section, __LOCATION__, name="ALMO_OPTIMIZER_TRUSTR", &
description="Controls the trust-region optimization of block-diagonal ALMOs. "// &
"See XALMO_OPTIMIZER_TRUSTR section for brief explanations.", &
n_keywords=10, n_subsections=0, repeats=.FALSE.)
optimizer_type = optimizer_trustr
CASE (optimizer_newton_pcg_solver)
CALL section_create(section, __LOCATION__, name="XALMO_NEWTON_PCG_SOLVER", &
description="Controls an iterative solver of the Newton-Raphson linear equation.", &
Expand Down Expand Up @@ -721,7 +734,7 @@ RECURSIVE SUBROUTINE create_optimizer_section(section, optimizer_id)
"how large the threshold is the maximum number of projected "// &
"eienvectors for a fragment equals to the number of occupied "// &
"orbitals of fragment's neighbors.", &
usage="PRECOND_FILTER_THRESHOLD 0.1", default_r_val=0.5_dp)
usage="PRECOND_FILTER_THRESHOLD 0.1", default_r_val=-1.0_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

Expand Down
1 change: 1 addition & 0 deletions tests/QS/regtest-almo-trustr/TEST_FILES
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
BN-cg.inp 11 1e-08 -101.701385613578893
FHchain-dogleg.inp 11 1e-06 -99.282920716902652
water-ionic-md-cauchy.inp 11 1e-08 -33.942336572891456
ice-block-diag-trustr.inp 11 1e-07 -137.557303714425842
#EOF
64 changes: 64 additions & 0 deletions tests/QS/regtest-almo-trustr/ice-block-diag-trustr.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
&GLOBAL
PROJECT hex-ice
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
POTENTIAL_FILE_NAME GTH_POTENTIALS
BASIS_SET_FILE_NAME GTH_BASIS_SETS
&QS
ALMO_SCF T
EPS_DEFAULT 1.0E-8 ! 1.0E-12
&END QS

&ALMO_SCF
EPS_FILTER 1.0E-8
ALMO_ALGORITHM TRUST_REGION
ALMO_SCF_GUESS ATOMIC
MO_OVERLAP_INV_ALG LS_HOTELLING

DELOCALIZE_METHOD NONE

&ALMO_OPTIMIZER_TRUSTR
ALGORITHM DOGLEG
MAX_ITER_OUTER_LOOP 50
EPS_ERROR 5.0E-4
INITIAL_TRUST_RADIUS 0.1
MAX_TRUST_RADIUS 2.0
&END ALMO_OPTIMIZER_TRUSTR

&END ALMO_SCF

&MGRID
CUTOFF 200 ! 320
NGRIDS 5
&END MGRID
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
@INCLUDE '../regtest-almo-1/ice.cell'
MULTIPLE_UNIT_CELL 1 1 1
&END CELL
&TOPOLOGY
MULTIPLE_UNIT_CELL 1 1 1
&END
&COORD
@INCLUDE '../regtest-almo-1/ice.xyz'
&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
&END SUBSYS
&END FORCE_EVAL

0 comments on commit 06b0715

Please sign in to comment.