higher orders (3rd-7th) of the Newton-Schulz-like sign iteration

cp2k · May 12, 2019 · 9d6c729 · 9d6c729
1 parent 61e0ea3
commit 9d6c729
Show file tree

Hide file tree

Showing 15 changed files with 830 additions and 33 deletions.
diff --git a/src/dm_ls_scf.F b/src/dm_ls_scf.F
@@ -637,8 +637,9 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env)
                   SELECT CASE (ls_scf_env%purification_method)
                   CASE (ls_scf_sign)
                      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, 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%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)
                   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), &
@@ -977,8 +978,9 @@ SUBROUTINE post_scf_mu_scan(ls_scf_env)
             IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
 
             CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, ls_scf_env%sign_method, &
-                                              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%sign_order, ls_scf_env%matrix_ks(ispin), &
+                                              ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, &
+                                              ls_scf_env%eps_filter)
          ENDDO
 
          t2 = m_walltime()

diff --git a/src/dm_ls_scf_create.F b/src/dm_ls_scf_create.F
@@ -30,7 +30,7 @@ MODULE dm_ls_scf_create
         ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, ls_s_inversion_none, &
         ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, &
         ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_pexsi, ls_scf_sign, &
-        ls_scf_tc2, ls_scf_trs4
+        ls_scf_sign_ns, ls_scf_sign_proot, ls_scf_tc2, ls_scf_trs4
    USE input_enumeration_types,         ONLY: enum_i2c,&
                                               enumeration_type
    USE input_keyword_types,             ONLY: keyword_get,&
@@ -210,6 +210,7 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
       CALL section_vals_val_get(ls_scf_section, "PERFORM_MU_SCAN", l_val=ls_scf_env%perform_mu_scan)
       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, "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)
@@ -383,6 +384,15 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
          SELECT CASE (ls_scf_env%purification_method)
          CASE (ls_scf_sign)
             WRITE (unit_nr, '(T2,A,T51,A30)') "Purification method", ADJUSTR("sign iteration")
+            SELECT CASE (ls_scf_env%sign_method)
+            CASE (ls_scf_sign_ns)
+               WRITE (unit_nr, '(T2,A,T61,A20)') "Sign method:", ADJUSTR("newton schulz")
+            CASE (ls_scf_sign_proot)
+               WRITE (unit_nr, '(T2,A,T61,A20)') "Sign method:", ADJUSTR("p-th root method")
+            CASE DEFAULT
+               CPABORT("Unkown sign method.")
+            END SELECT
+            WRITE (unit_nr, '(T2,A,T61,I20)') "Sign order:", ls_scf_env%sign_order
          CASE (ls_scf_tc2)
             CALL cite_reference(Niklasson2014)
             WRITE (unit_nr, '(T2,A,T51,A30)') "Purification method", ADJUSTR("Trace conserving 2nd order")

diff --git a/src/dm_ls_scf_methods.F b/src/dm_ls_scf_methods.F
@@ -365,6 +365,7 @@ END SUBROUTINE apply_matrix_preconditioner
 !> \param mu ...
 !> \param fixed_mu ...
 !> \param sign_method ...
+!> \param sign_order ...
 !> \param matrix_ks ...
 !> \param matrix_s ...
 !> \param matrix_s_inv ...
@@ -374,12 +375,13 @@ END SUBROUTINE apply_matrix_preconditioner
 !>       2010.10 created [Joost VandeVondele]
 !> \author Joost VandeVondele
 ! **************************************************************************************************
-   SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, matrix_ks, matrix_s, matrix_s_inv, nelectron, threshold)
+   SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, &
+                                  matrix_s, matrix_s_inv, nelectron, threshold)
 
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_p
       REAL(KIND=dp), INTENT(INOUT)                       :: mu
       LOGICAL                                            :: fixed_mu
-      INTEGER                                            :: sign_method
+      INTEGER                                            :: 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
@@ -414,7 +416,8 @@ SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, matrix_ks, m
             IF (ABS(mu_high-mu_low) < threshold) EXIT
          ENDIF
 
-         CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, matrix_ks, matrix_s, matrix_s_inv, threshold)
+         CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, &
+                                           matrix_ks, matrix_s, matrix_s_inv, threshold)
          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
 
@@ -445,6 +448,7 @@ END SUBROUTINE density_matrix_sign
 !> \param trace ...
 !> \param mu ...
 !> \param sign_method ...
+!> \param sign_order ...
 !> \param matrix_ks ...
 !> \param matrix_s ...
 !> \param matrix_s_inv ...
@@ -453,12 +457,13 @@ END SUBROUTINE density_matrix_sign
 !>       2010.10 created [Joost VandeVondele]
 !> \author Joost VandeVondele
 ! **************************************************************************************************
-   SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, matrix_ks, matrix_s, matrix_s_inv, threshold)
+   SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, matrix_ks, &
+                                           matrix_s, matrix_s_inv, threshold)
 
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_p
       REAL(KIND=dp), INTENT(OUT)                         :: trace
       REAL(KIND=dp), INTENT(INOUT)                       :: mu
-      INTEGER                                            :: sign_method
+      INTEGER                                            :: sign_method, sign_order
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_ks, matrix_s, matrix_s_inv
       REAL(KIND=dp), INTENT(IN)                          :: threshold
 
@@ -490,9 +495,9 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, matrix
       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)
+         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)
+         CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
       CASE DEFAULT
          CPABORT("Unkown sign method.")
       END SELECT

diff --git a/src/dm_ls_scf_types.F b/src/dm_ls_scf_types.F
@@ -102,6 +102,7 @@ MODULE dm_ls_scf_types
       INTEGER :: s_inversion_type
       INTEGER :: purification_method
       INTEGER :: sign_method
+      INTEGER :: sign_order
       INTEGER :: s_sqrt_method
       INTEGER :: s_sqrt_order
 

diff --git a/src/input_cp2k_ls.F b/src/input_cp2k_ls.F
@@ -215,6 +215,13 @@ SUBROUTINE create_ls_scf_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, name="SIGN_ORDER", &
+                          description="Order of the method used for the computation of the sign matrix.", &
+                          usage="SIGN_ORDER 2", &
+                          default_i_val=2)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, 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.)