Framework to calculate response forces for external energy expression…

…s based on KS orbitals (cp2k#3721)

src/CMakeLists.txt

@@ -134,6 +134,7 @@ list(
   ec_efield_local.F
   ec_environment.F
   ec_env_types.F
+  ec_external.F
   ec_methods.F
   ec_orth_solver.F
   ed_analysis.F
@@ -546,7 +547,6 @@ list(
   qs_elf_methods.F
   qs_energy.F
   qs_energy_init.F
-  qs_energy_matrix_w.F
   qs_energy_types.F
   qs_energy_utils.F
   qs_energy_window.F

src/ec_env_types.F

@@ -14,6 +14,8 @@
 MODULE ec_env_types
    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
    USE cp_dbcsr_operations,             ONLY: dbcsr_deallocate_matrix_set
+   USE cp_fm_types,                     ONLY: cp_fm_release,&
+                                              cp_fm_type
    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type,&
                                               ls_scf_release
    USE hfx_types,                       ONLY: hfx_release,&
@@ -66,6 +68,7 @@ MODULE ec_env_types
       ! debug
       LOGICAL                                          :: debug_forces = .FALSE.
       LOGICAL                                          :: debug_stress = .FALSE.
+      LOGICAL                                          :: debug_external = .FALSE.
       ! basis set
       CHARACTER(len=20)                                :: basis = ""
       LOGICAL                                          :: mao = .FALSE.
@@ -101,6 +104,9 @@ MODULE ec_env_types
                                                           matrix_w => Null()
       ! reduce basis
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER        :: mao_coef => Null()
+      ! external energy calclulation
+      TYPE(cp_fm_type), DIMENSION(:), POINTER          :: mo_occ => NULL()
+      TYPE(cp_fm_type), DIMENSION(:), POINTER          :: cpmos => NULL()
       ! CP equations
       TYPE(qs_p_env_type), POINTER                     :: p_env => Null()
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER        :: matrix_hz => Null(), matrix_wz => Null(), &
@@ -206,6 +212,9 @@ SUBROUTINE ec_env_release(ec_env)
          IF (ASSOCIATED(ec_env%xc_section_aux)) CALL section_vals_release(ec_env%xc_section_aux)
          IF (ASSOCIATED(ec_env%xc_section_primary)) CALL section_vals_release(ec_env%xc_section_primary)
 
+         CALL cp_fm_release(ec_env%mo_occ)
+         CALL cp_fm_release(ec_env%cpmos)
+
          DEALLOCATE (ec_env)
 
       END IF

src/ec_environment.F

@@ -28,12 +28,12 @@ MODULE ec_environment
    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
    USE ec_env_types,                    ONLY: energy_correction_type
    USE input_constants,                 ONLY: &
-        ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, &
-        ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, kg_cholesky, 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, xc_vdw_fun_nonloc, &
-        xc_vdw_fun_pairpot
+        ec_diagonalization, ec_functional_dc, ec_functional_ext, ec_functional_harris, &
+        ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, &
+        kg_cholesky, 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, &
+        xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
    USE input_cp2k_check,                ONLY: xc_functionals_expand
    USE input_section_types,             ONLY: section_get_ival,&
                                               section_vals_get,&
@@ -198,6 +198,8 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section)
                                    l_val=ec_env%debug_forces)
          CALL section_vals_val_get(ec_section, "DEBUG_STRESS", &
                                    l_val=ec_env%debug_stress)
+         CALL section_vals_val_get(ec_section, "DEBUG_EXTERNAL_METHOD", &
+                                   l_val=ec_env%debug_external)
          ! ADMM
          CALL section_vals_val_get(ec_section, "ADMM", l_val=ec_env%do_ec_admm)
 
@@ -281,13 +283,20 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section)
                           "DC-DFT: Correction and ground state need to use the same basis. "// &
                           "Checked by comparing basis set names only.")
          END IF
+         IF (ec_env%energy_functional == ec_functional_ext .AND. ec_env%basis_inconsistent) THEN
+            CALL cp_abort(__LOCATION__, &
+                          "Exteranl Energy: Correction and ground state need to use the same basis. "// &
+                          "Checked by comparing basis set names only.")
+         END IF
          !
          ! set functional
          SELECT CASE (ec_env%energy_functional)
          CASE (ec_functional_harris)
             ec_env%ec_name = "Harris"
          CASE (ec_functional_dc)
             ec_env%ec_name = "DC-DFT"
+         CASE (ec_functional_ext)
+            ec_env%ec_name = "External Energy"
          CASE DEFAULT
             CPABORT("unknown energy correction")
          END SELECT
@@ -502,6 +511,8 @@ SUBROUTINE ec_write_input(ec_env)
             WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL"
          CASE (ec_functional_dc)
             WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT"
+         CASE (ec_functional_ext)
+            WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "External"
          END SELECT
          WRITE (unit_nr, '()')