Skip to content

Commit

Permalink
LS: Add symmetric orthogonalization
Browse files Browse the repository at this point in the history 
Right now H is orthogonalized by multiplying with S^(-1), which does not
retain symmetry. Add an optional symmetric orthogonalization by
multiplying with S^(-1/2) from both sides.

This method can be activated by setting SIGN_SYMMETRIC in the input file
which defaults to .FALSE.

Co-authored-by: Robert Schade <robert.schade@uni-paderborn.de>
  • Loading branch information
@michaellass @robertschade @oschuett
2 people authored and oschuett committed Mar 12, 2020
1 parent c739195 commit c2469e9
Show file tree
Hide file tree
Showing 6 changed files with 85 additions and 27 deletions.
4 changes: 2 additions & 2 deletions src/dm_ls_scf.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -639,7 +639,7 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env)
CALL density_matrix_sign(ls_scf_env%matrix_p(ispin), ls_scf_env%mu_spin(ispin), ls_scf_env%fixed_mu, &
ls_scf_env%sign_method, ls_scf_env%sign_order, matrix_mixing_old(ispin), &
ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, nelectron_spin_real, &
ls_scf_env%eps_filter)
ls_scf_env%eps_filter, ls_scf_env%sign_symmetric)
CASE (ls_scf_tc2)
CALL density_matrix_tc2(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
Expand Down Expand Up @@ -980,7 +980,7 @@ SUBROUTINE post_scf_mu_scan(ls_scf_env)
CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, ls_scf_env%sign_method, &
ls_scf_env%sign_order, ls_scf_env%matrix_ks(ispin), &
ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, &
ls_scf_env%eps_filter)
ls_scf_env%eps_filter, ls_scf_env%sign_symmetric)
ENDDO

t2 = m_walltime()
Expand Down
2 changes: 2 additions & 0 deletions src/dm_ls_scf_create.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,6 +211,7 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
CALL section_vals_val_get(ls_scf_section, "PURIFICATION_METHOD", i_val=ls_scf_env%purification_method)
CALL section_vals_val_get(ls_scf_section, "SIGN_METHOD", i_val=ls_scf_env%sign_method)
CALL section_vals_val_get(ls_scf_section, "SIGN_ORDER", i_val=ls_scf_env%sign_order)
CALL section_vals_val_get(ls_scf_section, "SIGN_SYMMETRIC", l_val=ls_scf_env%sign_symmetric)
CALL section_vals_val_get(ls_scf_section, "DYNAMIC_THRESHOLD", l_val=ls_scf_env%dynamic_threshold)
CALL section_vals_val_get(ls_scf_section, "NON_MONOTONIC", l_val=ls_scf_env%non_monotonic)
CALL section_vals_val_get(ls_scf_section, "S_SQRT_METHOD", i_val=ls_scf_env%s_sqrt_method)
Expand Down Expand Up @@ -393,6 +394,7 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
CPABORT("Unkown sign method.")
END SELECT
WRITE (unit_nr, '(T2,A,T61,I20)') "Sign order:", ls_scf_env%sign_order
WRITE (unit_nr, '(T2,A,T61,L20)') "Symmetric sign calculation:", ls_scf_env%sign_symmetric
CASE (ls_scf_tc2)
CALL cite_reference(Niklasson2014)
WRITE (unit_nr, '(T2,A,T51,A30)') "Purification method", ADJUSTR("Trace conserving 2nd order")
Expand Down
96 changes: 72 additions & 24 deletions src/dm_ls_scf_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -372,12 +372,13 @@ END SUBROUTINE apply_matrix_preconditioner
!> \param matrix_s_inv ...
!> \param nelectron ...
!> \param threshold ...
!> \param sign_symmetric ...
!> \par History
!> 2010.10 created [Joost VandeVondele]
!> \author Joost VandeVondele
! **************************************************************************************************
SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, &
matrix_s, matrix_s_inv, nelectron, threshold)
matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric)

TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
REAL(KIND=dp), INTENT(INOUT) :: mu
Expand All @@ -386,6 +387,7 @@ SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order,
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s, matrix_s_inv
INTEGER, INTENT(IN) :: nelectron
REAL(KIND=dp), INTENT(IN) :: threshold
LOGICAL :: sign_symmetric

CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign', &
routineP = moduleN//':'//routineN
Expand Down Expand Up @@ -418,7 +420,8 @@ SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order,
ENDIF

CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, &
matrix_ks, matrix_s, matrix_s_inv, threshold)
matrix_ks, matrix_s, matrix_s_inv, threshold, &
sign_symmetric)
IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,I2,1X,F13.9,1X,F15.9)') &
"Density matrix: iter, mu, trace error: ", iter, mu, trace - nelectron

Expand Down Expand Up @@ -454,28 +457,31 @@ END SUBROUTINE density_matrix_sign
!> \param matrix_s ...
!> \param matrix_s_inv ...
!> \param threshold ...
!> \param sign_symmetric ...
!> \par History
!> 2010.10 created [Joost VandeVondele]
!> \author Joost VandeVondele
! **************************************************************************************************
SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, matrix_ks, &
matrix_s, matrix_s_inv, threshold)
matrix_s, matrix_s_inv, threshold, sign_symmetric)

TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
REAL(KIND=dp), INTENT(OUT) :: trace
REAL(KIND=dp), INTENT(INOUT) :: mu
INTEGER :: sign_method, sign_order
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s, matrix_s_inv
REAL(KIND=dp), INTENT(IN) :: threshold
LOGICAL :: sign_symmetric

CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign_fixed_mu', &
routineP = moduleN//':'//routineN

INTEGER :: handle, unit_nr
LOGICAL :: converged = .FALSE., symmetrize = .TRUE.
REAL(KIND=dp) :: frob_matrix
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_type) :: matrix_p_ud, matrix_sign, &
matrix_sinv_ks, matrix_tmp
TYPE(dbcsr_type) :: matrix_p_ud, matrix_s_sqrt, matrix_s_sqrt_inv, matrix_sign, &
matrix_sinv_ks, matrix_ssqrtinv_ks_ssqrtinv, matrix_ssqrtinv_ks_ssqrtinv2, matrix_tmp

CALL timeset(routineN, handle)

Expand All @@ -486,23 +492,56 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
unit_nr = -1
ENDIF

! get inv(S)*H-I*mu
CALL dbcsr_create(matrix_sinv_ks, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_inv, matrix_ks, &
0.0_dp, matrix_sinv_ks, filter_eps=threshold)
CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu)

! compute sign(inv(S)*H-I*mu)
CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
SELECT CASE (sign_method)
CASE (ls_scf_sign_ns)
CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE (ls_scf_sign_proot)
CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE DEFAULT
CPABORT("Unkown sign method.")
END SELECT
CALL dbcsr_release(matrix_sinv_ks)

IF (sign_symmetric) THEN
CALL dbcsr_create(matrix_s_sqrt, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(matrix_s_sqrt_inv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
CALL matrix_sqrt_Newton_Schulz(matrix_s_sqrt, matrix_s_sqrt_inv, matrix_s, threshold, sign_order, &
1.0E-10_dp, 100, symmetrize, converged) !FIXME
CALL dbcsr_release(matrix_s_sqrt)
IF (.NOT. converged) THEN
PRINT *, "sqrt inv did not converge"
STOP
ENDIF

CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)

SELECT CASE (sign_method)
CASE (ls_scf_sign_ns)
CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
CASE (ls_scf_sign_proot)
CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
CASE DEFAULT
CPABORT("Unkown sign method.")
END SELECT
CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)

ELSE ! .NOT. sign_symmetric
! get inv(S)*H-I*mu
CALL dbcsr_create(matrix_sinv_ks, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_inv, matrix_ks, &
0.0_dp, matrix_sinv_ks, filter_eps=threshold)
CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu)

! compute sign(inv(S)*H-I*mu)
SELECT CASE (sign_method)
CASE (ls_scf_sign_ns)
CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE (ls_scf_sign_proot)
CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE DEFAULT
CPABORT("Unkown sign method.")
END SELECT
CALL dbcsr_release(matrix_sinv_ks)
ENDIF

! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
Expand All @@ -523,9 +562,18 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
CALL dbcsr_release(matrix_tmp)
IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix

! get P=PS*inv(S)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_s_inv, &
0.0_dp, matrix_p, filter_eps=threshold)
IF (sign_symmetric) THEN
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
0.0_dp, matrix_p, filter_eps=threshold)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, matrix_s_sqrt_inv, &
0.0_dp, matrix_p, filter_eps=threshold)
CALL dbcsr_release(matrix_s_sqrt_inv)
ELSE

! get P=PS*inv(S)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_s_inv, &
0.0_dp, matrix_p, filter_eps=threshold)
ENDIF
CALL dbcsr_release(matrix_p_ud)

CALL timestop(handle)
Expand Down
1 change: 1 addition & 0 deletions src/dm_ls_scf_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,7 @@ MODULE dm_ls_scf_types
INTEGER :: purification_method
INTEGER :: sign_method
INTEGER :: sign_order
LOGICAL :: sign_symmetric
INTEGER :: s_sqrt_method
INTEGER :: s_sqrt_order

Expand Down
3 changes: 2 additions & 1 deletion src/energy_corrections.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1542,7 +1542,8 @@ SUBROUTINE ec_ls_solver(qs_env, ls_method, matrix_ks, matrix_s, matrix_p, matrix
CASE (ec_matrix_sign)
CALL density_matrix_sign(matrix_p(ispin, 1)%matrix, mu_spin(ispin), fixed_mu, &
sign_method, sign_order, matrix_ks(ispin, 1)%matrix, &
matrix_s(1, 1)%matrix, matrix_s_inv, nelectron, eps_filter)
matrix_s(1, 1)%matrix, matrix_s_inv, nelectron, eps_filter, &
sign_symmetric=.FALSE.)
CASE (ec_matrix_trs4)
CALL density_matrix_trs4(matrix_p(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix, &
matrix_s_sqrt_inv, nelectron, eps_filter, &
Expand Down
6 changes: 6 additions & 0 deletions src/input_cp2k_ls.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,12 @@ SUBROUTINE create_ls_scf_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="SIGN_SYMMETRIC", &
description="Use symmetric orthogonalization when generating the input for the sign function.", &
usage="SIGN_SYMMETRIC .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="DYNAMIC_THRESHOLD", &
description="Should the threshold for the purification be chosen dynamically", &
usage="DYNAMIC_THRESHOLD .TRUE.", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
Expand Down

0 comments on commit c2469e9

Please sign in to comment.