Optionally switch off ADMM XC correction in TDDFPT kernel (#2205)

* GAPW_XC bug fix for linear response

* GAPW/GAPW_XC linear response refactoring

* Prettify and update regtests

* Remove filter and adjust regtest value

* Pretty

* Linear response for ADMM/GAPW

* TDDFT/ADMM/GAPW excitation energies. More regtests. Some refactoring.

* Merge with upstream

* Optionally switch off ADMM XC correction in TDDFPT kernel

src/cp_control_types.F

@@ -475,6 +475,8 @@ MODULE cp_control_types
       INTEGER                              :: auto_basis_p_lri_aux = 1
       !> use symmetric definition of ADMM Kernel correction
       LOGICAL                              :: admm_symm
+      !> Use/Ignore possible ADMM Kernel XC correction
+      LOGICAL                              :: admm_xc_correction
       !
       ! DIPOLE_MOMENTS subsection
       !

src/cp_control_utils.F

@@ -1282,6 +1282,7 @@ SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
       CALL section_vals_val_get(t_section, "RKS_TRIPLETS", l_val=t_control%rks_triplets)
       CALL section_vals_val_get(t_section, "DO_LRIGPW", l_val=t_control%do_lrigpw)
       CALL section_vals_val_get(t_section, "ADMM_KERNEL_CORRECTION_SYMMETRIC", l_val=t_control%admm_symm)
+      CALL section_vals_val_get(t_section, "ADMM_KERNEL_XC_CORRECTION", l_val=t_control%admm_xc_correction)
 
       ! read automatically generated auxiliary basis for LRI
       CALL section_vals_val_get(t_section, "AUTO_BASIS", n_rep_val=nrep)

src/input_cp2k_properties_dft.F

@@ -1386,6 +1386,14 @@ SUBROUTINE create_tddfpt2_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="ADMM_KERNEL_XC_CORRECTION", &
+                          description="Use/Ignore ADMM correction xc functional for TD kernel."// &
+                          "XC correction functional is defined in ground state XC section.", &
+                          n_var=1, type_of_var=logical_t, &
+                          default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="ADMM_KERNEL_CORRECTION_SYMMETRIC", &
                           description="ADMM correction functional in kernel is applied symmetrically."// &
                           "Original implementation is using a non-symmetric formula.", &

src/qs_tddfpt2_eigensolver.F

@@ -364,7 +364,8 @@ SUBROUTINE tddfpt_compute_Aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfp
             ! full TDDFPT kernel
             CALL fhxc_kernel(Aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, &
                              kernel_env%full_kernel, kernel_env_admm_aux, sub_env, work_matrices, &
-                             tddfpt_control%admm_symm, do_lri_response)
+                             tddfpt_control%admm_symm, tddfpt_control%admm_xc_correction, &
+                             do_lri_response)
          ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
             ! sTDA kernel
             CALL stda_kernel(Aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, &

src/qs_tddfpt2_fhxc_forces.F

@@ -406,7 +406,7 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
          END IF
          CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
 
-         IF (do_admm) THEN
+         IF (do_admm .AND. tddfpt_control%admm_xc_correction) THEN
             IF (dft_control%admm_control%aux_exch_func == do_admm_aux_exch_func_none) THEN
                ! nothing to do
             ELSE

src/qs_tddfpt2_methods.F

@@ -269,10 +269,17 @@ SUBROUTINE tddfpt(qs_env, calc_forces)
       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
          IF (do_admm) THEN
             ! Full kernel with ADMM
-            CALL create_kernel_env(kernel_env=full_kernel_env, &
-                                   rho_struct_sub=work_matrices%rho_orb_struct_sub, &
-                                   xc_section=admm_env%xc_section_primary, &
-                                   is_rks_triplets=tddfpt_control%rks_triplets, sub_env=sub_env)
+            IF (tddfpt_control%admm_xc_correction) THEN
+               CALL create_kernel_env(kernel_env=full_kernel_env, &
+                                      rho_struct_sub=work_matrices%rho_orb_struct_sub, &
+                                      xc_section=admm_env%xc_section_primary, &
+                                      is_rks_triplets=tddfpt_control%rks_triplets, sub_env=sub_env)
+            ELSE
+               CALL create_kernel_env(kernel_env=full_kernel_env, &
+                                      rho_struct_sub=work_matrices%rho_orb_struct_sub, &
+                                      xc_section=xc_section, &
+                                      is_rks_triplets=tddfpt_control%rks_triplets, sub_env=sub_env)
+            END IF
 
             CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, &
                                                   rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, &
@@ -869,12 +876,15 @@ SUBROUTINE kernel_info(log_unit, dft_control, tddfpt_control, xc_section)
             CALL xc_write(log_unit, xc_section, lsd)
             IF (tddfpt_control%do_hfx) THEN
                IF (tddfpt_control%do_admm) THEN
-                  WRITE (log_unit, "(T2,A,T62,A19)") "KERNEL|", "ADMM exact exchange"
+                  WRITE (log_unit, "(T2,A,T62,A19)") "KERNEL|", "ADMM Exact Exchange"
+                  IF (tddfpt_control%admm_xc_correction) THEN
+                     WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Apply ADMM Kernel XC Correction"
+                  END IF
                   IF (tddfpt_control%admm_symm) THEN
-                     WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "symmetric ADMM kernel"
+                     WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Symmetric ADMM Kernel"
                   END IF
                ELSE
-                  WRITE (log_unit, "(T2,A,T67,A14)") "KERNEL|", "exact exchange"
+                  WRITE (log_unit, "(T2,A,T67,A14)") "KERNEL|", "Exact Exchange"
                END IF
             END IF
             IF (tddfpt_control%do_lrigpw) THEN

tests/QS/regtest-debug-1/TEST_FILES

@@ -10,5 +10,4 @@ h2o_periodic.inp                                      87      1e-05
 h2o_gga.inp                                           87      1e-05             0.165143521103E+02
 h2o_gapw.inp                                          87      1e-05             0.170288288429E+02
 h2o_gapw_xc.inp                                       87      1e-05             0.170142638188E+02
-h2o_admm_gapw.inp                                     87      1e-05             0.165349668089E+02
 #EOF

tests/QS/regtest-debug-2/TEST_FILES

@@ -7,4 +7,5 @@ h2o_lri.inp                                           87      1e-05
 h2o_pade_fd.inp                                       87      1e-05             0.168333363169E+02
 h2o_hfx.inp                                           87      1e-05             0.155040875037E+02
 h2o_hfx_admm.inp                                      87      1e-05             0.153251090463E+02
+h2o_admm_gapw.inp                                     87      1e-05             0.165349668089E+02
 #EOF

tests/QS/regtest-debug-4/TEST_FILES

@@ -16,4 +16,5 @@ h2o_dip09.inp                                          86      1e-09
 h2o_dip10.inp                                          86      1e-09             0.150644266924E+00
 h2o_dip11.inp                                          86      1e-09             0.804659694314E-01
 h2o_dip11b.inp                                         86      1e-09             0.104678363910E+00
+h2o_dip12.inp                                          86      1e-09            -0.129459747287E+00
 #EOF

tests/QS/regtest-debug-4/h2o_dip12.inp

@@ -0,0 +1,96 @@
+&FORCE_EVAL
+  METHOD Quickstep
+  &PROPERTIES
+    &TDDFPT
+       KERNEL FULL
+       ADMM_KERNEL_CORRECTION_SYMMETRIC   T
+       ADMM_KERNEL_XC_CORRECTION          F
+       NSTATES      1
+       MAX_ITER    50
+       CONVERGENCE [eV] 1.0e-7
+    &END TDDFPT
+  &END PROPERTIES
+  &DFT
+    &QS
+      METHOD GPW
+    &END QS
+    BASIS_SET_FILE_NAME BASIS_SET
+    BASIS_SET_FILE_NAME BASIS_ADMM
+    &AUXILIARY_DENSITY_MATRIX_METHOD
+      ADMM_PURIFICATION_METHOD NONE
+      EXCH_CORRECTION_FUNC PBEX
+      EXCH_SCALING_MODEL NONE
+      METHOD BASIS_PROJECTION
+    &END 
+    &EFIELD
+    &END
+    &EXCITED_STATES T
+       STATE 1
+    &END EXCITED_STATES
+    &SCF
+      SCF_GUESS ATOMIC 
+      &OT
+         PRECONDITIONER FULL_SINGLE_INVERSE
+         MINIMIZER DIIS
+      &END
+      &OUTER_SCF
+         MAX_SCF  20
+         EPS_SCF 1.0E-7
+      &END
+      MAX_SCF  50
+      EPS_SCF 1.0E-7
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE0
+      &END XC_FUNCTIONAL
+    &END XC
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+    &POISSON
+       PERIODIC NONE
+       POISSON_SOLVER MT
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 4.0 4.0 4.0
+      PERIODIC NONE
+    &END
+    &COORD
+    O   0.000000    0.000000   -0.065587
+    H   0.000000   -0.757136    0.520545
+    H   0.000000    0.757136    0.520545
+    &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END
+    &KIND H
+      BASIS_SET DZV-GTH-PADE
+      BASIS_SET AUX_FIT FIT3
+      POTENTIAL GTH-PADE-q1
+    &END KIND
+    &KIND O
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT FIT3
+      POTENTIAL GTH-PADE-q6
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT dipole
+  RUN_TYPE DEBUG       
+&END GLOBAL
+&DEBUG
+   DEBUG_FORCES .FALSE.
+   DEBUG_STRESS_TENSOR .FALSE.
+   DEBUG_DIPOLE .TRUE.
+   DEBUG_POLARIZABILITY .FALSE.
+   DE 0.0002
+&END
+

tests/QS/regtest-tddfpt-admm/TEST_FILES

@@ -0,0 +1,11 @@
+# runs are executed in the same order as in this file
+# the second field tells which test should be run in order to compare with the last available output
+# e.g. 0 means do not compare anything, running is enough
+#      1 compares the last total energy in the file
+#      for details see cp2k/tools/do_regtest
+#
+h2o_a01.inp                                                1    5.0E-11             -17.21071638924079
+h2o_a02.inp                                                1    5.0E-11             -17.20961421163643
+h2o_a03.inp                                                1    5.0E-11             -17.20804254307528
+h2o_a04.inp                                                1    5.0E-11             -17.20949405791527
+#EOF

tests/QS/regtest-tddfpt-admm/h2o_a01.inp

@@ -0,0 +1,107 @@
+&FORCE_EVAL
+  METHOD Quickstep
+  &PROPERTIES
+    &TDDFPT
+       KERNEL FULL
+       ADMM_KERNEL_CORRECTION_SYMMETRIC   F
+       ADMM_KERNEL_XC_CORRECTION          F
+       &XC
+         &XC_FUNCTIONAL PBE
+         &END XC_FUNCTIONAL
+       &END XC
+       NSTATES      2
+       MAX_ITER    50
+       CONVERGENCE [eV] 1.0e-7
+       RKS_TRIPLETS F
+    &END TDDFPT
+  &END PROPERTIES
+  &DFT
+    &QS
+      METHOD GPW
+    &END QS
+    BASIS_SET_FILE_NAME BASIS_SET
+    BASIS_SET_FILE_NAME BASIS_ADMM
+    &AUXILIARY_DENSITY_MATRIX_METHOD
+      ADMM_PURIFICATION_METHOD NONE
+      EXCH_CORRECTION_FUNC PBEX
+      EXCH_SCALING_MODEL NONE
+      METHOD BASIS_PROJECTION
+    &END
+    &EXCITED_STATES T
+       STATE 1
+    &END EXCITED_STATES
+    &SCF
+      SCF_GUESS ATOMIC 
+      &OT
+         PRECONDITIONER FULL_SINGLE_INVERSE
+         MINIMIZER DIIS
+      &END
+      &OUTER_SCF
+         MAX_SCF  20
+         EPS_SCF 1.0E-7
+      &END
+      MAX_SCF  50
+      EPS_SCF 1.0E-7
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE0
+      &END XC_FUNCTIONAL
+    &END XC
+    &MGRID
+       CUTOFF  300
+       REL_CUTOFF  60
+    &END
+    &POISSON
+       PERIODIC NONE
+       POISSON_SOLVER MT
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 4.0 4.0 4.0
+      PERIODIC NONE
+    &END
+    &COORD
+    O   0.000000    0.000000    0.000000
+    H   0.000000   -0.757136    0.580545
+    H   0.000000    0.757136    0.580545
+    &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END
+    &KIND H
+      BASIS_SET DZV-GTH-PADE
+      BASIS_SET AUX_FIT FIT3
+      POTENTIAL GTH-PADE-q1
+    &END KIND
+    &KIND O
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT FIT3
+      POTENTIAL GTH-PADE-q6
+    &END KIND
+  &END SUBSYS
+  &PRINT
+     &FORCES
+     &END
+  &END
+&END FORCE_EVAL
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT ftest
+  RUN_TYPE GEO_OPT
+&END GLOBAL
+&MOTION
+  &GEO_OPT
+     MAX_ITER 1
+  &END
+&END
+&DEBUG
+   DEBUG_FORCES .TRUE.
+   DEBUG_STRESS_TENSOR .FALSE.
+   DEBUG_DIPOLE .FALSE.
+   DEBUG_POLARIZABILITY .FALSE.
+   STOP_ON_MISMATCH F
+   DE 0.0005
+&END
+