Keyword to enable explicit switching off of sTDA exchange (#1805)

src/cp_control_types.F

@@ -420,6 +420,7 @@ MODULE cp_control_types
 ! **************************************************************************************************
    TYPE stda_control_type
       LOGICAL                                :: do_ewald
+      LOGICAL                                :: do_exchange
       REAL(KIND=dp)                          :: hfx_fraction
       REAL(KIND=dp)                          :: eps_td_filter
    END TYPE stda_control_type

src/cp_control_utils.F

@@ -43,7 +43,7 @@ MODULE cp_control_utils
         real_time_propagation, sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, &
         sccs_derivative_cd7, sccs_derivative_fft, sccs_fattebert_gygi, sic_ad, sic_eo, &
         sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, slater, &
-        tddfpt_excitations, use_mom_ref_user, xas_dip_len
+        tddfpt_excitations, tddfpt_kernel_stda, use_mom_ref_user, xas_dip_len
    USE input_cp2k_check,                ONLY: xc_functionals_expand
    USE input_cp2k_dft,                  ONLY: create_dft_section
    USE input_enumeration_types,         ONLY: enum_i2c,&
@@ -1243,7 +1243,9 @@ SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
 
       INTEGER                                            :: handle
       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: inds
-      LOGICAL                                            :: explicit, multigrid_set
+      LOGICAL                                            :: do_ewald, do_exchange, expl, explicit, &
+                                                            multigrid_set
+      REAL(KIND=dp)                                      :: filter, hfx
       TYPE(section_vals_type), POINTER                   :: dipole_section, mgrid_section, &
                                                             stda_section, xc_func, xc_section
 
@@ -1357,16 +1359,23 @@ SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
 
       ! sTDA subsection
       stda_section => section_vals_get_subs_vals(t_section, "STDA")
-      CALL section_vals_get(stda_section, explicit=explicit)
-      IF (explicit) THEN
-         CALL section_vals_val_get(stda_section, "HFX_FRACTION", r_val=t_control%stda_control%hfx_fraction)
-         CALL section_vals_val_get(stda_section, "EPS_TD_FILTER", r_val=t_control%stda_control%eps_td_filter)
-         CALL section_vals_val_get(stda_section, "DO_EWALD", l_val=t_control%stda_control%do_ewald)
-      ELSE
+      IF (t_control%kernel == tddfpt_kernel_stda) THEN
          t_control%stda_control%hfx_fraction = 0.0_dp
+         t_control%stda_control%do_exchange = .TRUE.
          t_control%stda_control%eps_td_filter = 1.e-10_dp
          t_control%stda_control%do_ewald = .FALSE.
       END IF
+      CALL section_vals_get(stda_section, explicit=explicit)
+      IF (explicit) THEN
+         CALL section_vals_val_get(stda_section, "HFX_FRACTION", r_val=hfx, explicit=expl)
+         IF (expl) t_control%stda_control%hfx_fraction = hfx
+         CALL section_vals_val_get(stda_section, "EPS_TD_FILTER", r_val=filter, explicit=expl)
+         IF (expl) t_control%stda_control%eps_td_filter = filter
+         CALL section_vals_val_get(stda_section, "DO_EWALD", l_val=do_ewald, explicit=expl)
+         IF (expl) t_control%stda_control%do_ewald = do_ewald
+         CALL section_vals_val_get(stda_section, "DO_EXCHANGE", l_val=do_exchange, explicit=expl)
+         IF (expl) t_control%stda_control%do_exchange = do_exchange
+      END IF
 
       CALL timestop(handle)
    END SUBROUTINE read_tddfpt2_control

src/input_cp2k_properties_dft.F

@@ -1598,6 +1598,14 @@ SUBROUTINE create_stda_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      ! even if scaling parameter for exchange FRACTION (see above) is zero, the semi-empirical electron repulsion
+      ! operator for exchange is not, so that a keyword is required to switch off sTDA exchange (if wanted)
+      CALL keyword_create(keyword, __LOCATION__, name="DO_EXCHANGE", &
+                          description="Explicitly including or switching off sTDA exchange", &
+                          usage="DO_EXCHANGE", 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="DO_EWALD", &
                           description="Use Ewald type method for Coulomb interaction", &
                           usage="DO_EWALD", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)

src/qs_tddfpt2_fhxc_forces.F

@@ -1231,7 +1231,7 @@ SUBROUTINE stda_force(qs_env, ex_env, gs_mos, stda_env, sub_env, work, debug_for
          !
          ! *** Exchange contribution
          !
-         IF (stda_env%hfx_fraction > 0.0_dp) THEN
+         IF (stda_env%do_exchange) THEN
             !
             IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
             !
@@ -1266,12 +1266,26 @@ SUBROUTINE stda_force(qs_env, ex_env, gs_mos, stda_env, sub_env, work, debug_for
                eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
                       stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
                rbeta = dr**bp
-               gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
+               IF (hfx > 0.0_dp) THEN
+                  gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
+               ELSE
+                  IF (dr < 1.0e-6_dp) THEN
+                     gabr = 0.0_dp
+                  ELSE
+                     gabr = 1._dp/dr
+                  END IF
+               END IF
+               !      gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
                ! forces
                IF (dr > 1.0e-6_dp) THEN
-                  dgabr = -(1._dp/bp)*(1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp + 1._dp)
-                  dgabr = bp*rbeta/dr**2*dgabr
-                  dgabr = SUM(pblock**2)*dgabr
+                  IF (hfx > 0.0_dp) THEN
+                     dgabr = -(1._dp/bp)*(1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp + 1._dp)
+                     dgabr = bp*rbeta/dr**2*dgabr
+                     dgabr = SUM(pblock**2)*dgabr
+                  ELSE
+                     dgabr = -1.0_dp/dr**3
+                     dgabr = SUM(pblock**2)*dgabr
+                  END IF
                   atom_i = atom_of_kind(iatom)
                   atom_j = atom_of_kind(jatom)
                   DO i = 1, 3

src/qs_tddfpt2_stda_types.F

@@ -39,8 +39,9 @@ MODULE qs_tddfpt2_stda_types
    END TYPE
 
    TYPE stda_env_type
-      !fraction of non-local HF exchange !! REDUNDANT also in stda_control
+      !fraction of non-local HF exchange
       REAL(KIND=dp)                                :: hfx_fraction
+      LOGICAL                                      :: do_exchange
       ! empirical parameters
       REAL(KIND=dp)                                :: alpha_param
       REAL(KIND=dp)                                :: beta_param
@@ -140,8 +141,15 @@ SUBROUTINE stda_init_param(qs_env, stda_kernel, stda_control)
                               z=stda_kernel%kind_param_set(ikind)%kind_param%z)
       END DO
 
-      ! HFx_fraction
-      stda_kernel%hfx_fraction = stda_control%hfx_fraction
+      IF (stda_control%do_exchange) THEN ! option to switch off exchange
+         ! HFx_fraction
+         stda_kernel%do_exchange = .TRUE.
+         stda_kernel%hfx_fraction = stda_control%hfx_fraction
+      ELSE
+         stda_kernel%do_exchange = .FALSE.
+         stda_kernel%hfx_fraction = 0.0_dp
+      END IF
+
       ! alpha and beta parameter
       stda_kernel%alpha_param = alpha(1) + stda_kernel%hfx_fraction*alpha(2)
       stda_kernel%beta_param = beta(1) + stda_kernel%hfx_fraction*beta(2)
@@ -179,9 +187,9 @@ SUBROUTINE allocate_stda_env(qs_env, stda_kernel, n_ao, nactive)
       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
       TYPE(stda_kind_type), POINTER                      :: kind_param
 
-      stda_kernel%hfx_fraction = 0.0d0
-      stda_kernel%alpha_param = 0.0d0
-      stda_kernel%beta_param = 0.0d0
+      stda_kernel%hfx_fraction = 0.0_dp
+      stda_kernel%alpha_param = 0.0_dp
+      stda_kernel%beta_param = 0.0_dp
       stda_kernel%nactive = 0
       stda_kernel%nactive(1:2) = nactive(1:2)
       stda_kernel%n_ao = n_ao

src/qs_tddfpt2_stda_utils.F

@@ -516,7 +516,7 @@ SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: first_sgf, kind_of, last_sgf
       INTEGER, DIMENSION(2)                              :: nactive, nlim
       LOGICAL                                            :: calculate_forces, do_coulomb, do_ewald, &
-                                                            use_virial
+                                                            do_exchange, use_virial
       REAL(KIND=dp)                                      :: alpha, bp, dr, eta, gabr, hfx, rbeta, &
                                                             spinfac
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: tcharge, tv
@@ -554,6 +554,7 @@ SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
          do_coulomb = .TRUE.
       END IF
       do_ewald = stda_control%do_ewald
+      do_exchange = stda_control%do_exchange
 
       para_env => sub_env%para_env
 
@@ -673,7 +674,8 @@ SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
          !
          ! *** Exchange contribution
          !
-         IF (stda_env%hfx_fraction > 0.0_dp) THEN
+         IF (do_exchange) THEN ! option to explicitly switch off exchange
+            ! (exchange contributes also if hfx_fraction=0)
             ALLOCATE (kind_of(natom))
             CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
             CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
@@ -702,7 +704,15 @@ SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
                eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
                       stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
                rbeta = dr**bp
-               gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
+               IF (hfx > 0.0_dp) THEN
+                  gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
+               ELSE
+                  IF (dr < 1.e-6) THEN
+                     gabr = 0.0_dp
+                  ELSE
+                     gabr = 1._dp/dr
+                  END IF
+               END IF
                pblock = gabr*pblock
             END DO
             CALL dbcsr_iterator_stop(iter)

tests/QS/regtest-debug-5/h2o_dip20.inp

@@ -5,6 +5,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-debug-5/h2o_dip21.inp

@@ -5,6 +5,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-debug-5/h2o_dip22.inp

@@ -5,6 +5,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-debug-5/h2o_dip25.inp

@@ -6,6 +6,7 @@
        RKS_TRIPLETS T
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-stda-force-3/h2o_f20.inp

@@ -5,6 +5,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.00
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-stda-force/h2o_f04.inp

@@ -5,6 +5,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.00
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-stda/CH2O_s1.inp

@@ -10,6 +10,7 @@
        KERNEL sTDA
        &sTDA 
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      5
        MAX_ITER    50

tests/QS/regtest-tddfpt-stda/CH2O_tddfpt_stda-pbe-s.inp

@@ -11,6 +11,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      1
        MAX_ITER    50

tests/QS/regtest-tddfpt-stda/CH2O_tddfpt_stda-pbe-s_doexchange.inp

@@ -0,0 +1,70 @@
+&GLOBAL
+  PROJECT CH2O_tddfpt_stda-s-1
+  RUN_TYPE ENERGY
+  PRINT_LEVEL LOW
+  PREFERRED_DIAG_LIBRARY SL
+&END GLOBAL
+&FORCE_EVAL
+  METHOD Quickstep
+  &PROPERTIES
+    &TDDFPT
+       KERNEL sTDA
+       &sTDA
+          DO_EXCHANGE T
+          FRACTION 0.0
+       &END sTDA
+       NSTATES      1
+       MAX_ITER    50
+       CONVERGENCE [eV] 1.0e-5
+    &END TDDFPT
+  &END PROPERTIES
+  &DFT
+    BASIS_SET_FILE_NAME GTH_BASIS_SETS
+    POTENTIAL_FILE_NAME POTENTIAL
+    &MGRID
+      CUTOFF 100
+    &END MGRID
+    &QS
+    &END QS
+    &POISSON
+     POISSON_SOLVER WAVELET
+     PERIODIC NONE
+    &END POISSON
+    &SCF
+      MAX_SCF 20
+      SCF_GUESS atomic
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE
+      &END XC_FUNCTIONAL
+    &END XC
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 8.0 8.0 8.0
+      PERIODIC NONE
+    &END CELL
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END TOPOLOGY
+    &COORD
+  O         0.094933   -0.000368    0.895642
+  C        -0.031077   -0.000121   -0.307326
+  H        -0.090437    0.947608   -0.895642
+  H        -0.094933   -0.947608   -0.895562
+    &END COORD
+    &KIND O
+      BASIS_SET aug-TZV2P-GTH
+      POTENTIAL GTH-PBE-q6
+    &END KIND
+    &KIND C
+      BASIS_SET aug-TZV2P-GTH
+      POTENTIAL GTH-PBE-q4
+    &END KIND
+    &KIND H
+      BASIS_SET aug-TZV2P-GTH
+      POTENTIAL GTH-PBE-q1
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL

tests/QS/regtest-tddfpt-stda/CH2O_tddfpt_stda-s-1.inp

@@ -10,6 +10,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      3
        MAX_ITER    500

tests/QS/regtest-tddfpt-stda/CH2O_tddfpt_stda-xtb-s.inp

@@ -11,6 +11,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      1
        MAX_ITER    500

tests/QS/regtest-tddfpt-stda/H2O_tddfpt_stda-pbe-s.inp

@@ -10,6 +10,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      3
        MAX_ITER    50

tests/QS/regtest-tddfpt-stda/H2O_tddfpt_stda-pbe-t.inp

@@ -10,6 +10,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      3
        RKS_TRIPLETS T

tests/QS/regtest-tddfpt-stda/NO_tddfpt_stda-s-1.inp

@@ -10,6 +10,7 @@
        KERNEL sTDA
        &sTDA
           FRACTION 0.0
+          DO_EXCHANGE F
        &END sTDA
        NSTATES      3
        MAX_ITER    50

tests/QS/regtest-tddfpt-stda/TEST_FILES

@@ -6,6 +6,7 @@
 CH2O_tddfpt_stda-s-1.inp                                   37    1.0E-05               0.349917E+00
 CH2Oplus_tddfpt_stda-lsd.inp                               37    1.0E-05               0.212156E+00
 CH2O_tddfpt_stda-pbe-s.inp                                 37    1.0E-05               0.125452E+00
+CH2O_tddfpt_stda-pbe-s_doexchange.inp                      37    1.0E-05               0.155105E+00
 CH2O_tddfpt_stda-xtb-s.inp                                 37    1.0E-05               0.150482E+00
 CH2O_stda-xtb-s.inp                                        37    1.0E-05               0.256889E-01
 H2O_tddfpt_stda-pbe-s.inp                                  37    1.0E-05               0.599658E+00