GAPW: towards response and TDDFT forces (#2345)

* TDDFT RHS of CPHF equation using GAPW Refectoring of ks_ref_potential and exstate_type * TDDFT GAPW forces: RHS CPHF and direct forces * Further GAPW/TDDFT Gradient terms ADMM Slater X correction * Update to new pointer scheme * Remove duplicated code and debug statements * Pretty * update data * Deallocate PW structure from qs_vxc in ks_reference * Get rid of memory leaks * rebase
cp2k · Oct 13, 2022 · cea20fb · cea20fb
1 parent 16f0f40
commit cea20fb
Show file tree

Hide file tree

Showing 30 changed files with 2,062 additions and 407 deletions.
diff --git a/data/POTENTIAL_UZH b/data/POTENTIAL_UZH
@@ -1,4 +1,164 @@
 ################################################################################
+#
+# All-electron (full nuclear potential) definitions in the GTH-format:
+#
+# Element symbol  Name of the potential  Alias names
+# n_elec(s)  n_elec(p)  n_elec(d)  ...
+# r_loc   0
+#
+# n_elec   : Number of electrons for each angular momentum quantum number
+#            (electronic configuration -> s p d ...)
+# r_loc    : Radius for the local part defined by the Gaussian function
+#            exponent alpha_erf
+#
+################################################################################
+#
+H ALLELECTRON ALL
+    1    0    0
+     0.20000000    0
+#
+He ALLELECTRON ALL
+    2    0    0
+     0.20000000    0
+#
+Li ALLELECTRON ALL
+    3    0    0
+     0.40000000    0
+#
+Be ALLELECTRON ALL
+    4    0    0
+     0.32500000    0
+#
+B ALLELECTRON ALL
+    4    1    0
+     0.43392956    0
+#
+C ALLELECTRON ALL
+    4    2    0
+     0.34883045    0
+#
+N ALLELECTRON ALL
+    4    3    0
+     0.28917923    0
+#
+O ALLELECTRON ALL
+    4    4    0
+     0.24762086    0
+#
+F ALLELECTRON ALL
+    4    5    0
+     0.21852465    0
+#
+Ne ALLELECTRON ALL
+    4    6    0
+     0.19000000    0
+#
+Na ALLELECTRON ALL
+    5    6    0
+     0.24631780    0
+#
+Mg ALLELECTRON ALL
+    6    6    0
+     0.21094954    0
+#
+Al ALLELECTRON ALL
+    6    7    0
+     0.45000000    0
+#
+Si ALLELECTRON ALL
+    6    8    0
+     0.44000000    0
+#
+P ALLELECTRON ALL
+    6    9    0
+     0.43000000    0
+#
+S ALLELECTRON ALL
+    6   10    0
+     0.42000000    0
+#
+Cl ALLELECTRON ALL
+    6   11    0
+     0.41000000    0
+#
+Ar ALLELECTRON ALL
+    6   12    0
+     0.40000000    0
+#
+K ALLELECTRON ALL
+    7   12    0
+     0.40000000    0
+#
+Ca ALLELECTRON ALL
+    8   12    0
+     0.39000000    0
+#
+Sc ALLELECTRON ALL
+    8   12    1
+     0.38500000    0
+#
+Ti ALLELECTRON ALL
+    8   12    2
+     0.3800000    0
+#
+V  ALLELECTRON ALL
+    8   12    3
+     0.37500000    0
+#
+Cr ALLELECTRON ALL
+    8   12    4
+     0.37000000    0
+#
+Mn ALLELECTRON ALL
+    8   12    5
+     0.36500000    0
+#
+Fe  ALLELECTRON ALL
+    8   12    6
+     0.36000000    0
+#
+Co  ALLELECTRON ALL
+    8   12    7
+     0.35500000    0
+#
+Ni  ALLELECTRON ALL
+    8   12    8
+     0.35000000    0
+#
+Cu ALLELECTRON ALL
+    8   12    9
+     0.53000000    0
+#
+Zn ALLELECTRON ALL
+    8   12   10
+     0.34000000    0
+#
+Ga ALLELECTRON ALL
+    8   13   10
+     0.49000000    0
+#
+Ge ALLELECTRON ALL
+    8   14   10
+     0.54000000    0
+#
+As ALLELECTRON ALL
+    8   15   10
+     0.52000000    0
+#
+Se ALLELECTRON ALL
+    8   16   10
+     0.51000000    0
+#
+Br ALLELECTRON ALL
+    8   17   10
+     0.50000000    0
+#
+Kr ALLELECTRON ALL
+    8   18   10
+     0.50000000    0
+#
+################################################################################
+#
 #            GTH POTENTIALS OPTIMIZED USING THE UZH PROTOCOL
 ################################################################################
 #

diff --git a/src/ec_env_types.F b/src/ec_env_types.F
@@ -101,7 +101,7 @@ MODULE ec_env_types
       ! Harris (rhoout), and response density (rhoz) on grid
       TYPE(pw_p_type), DIMENSION(:), POINTER           :: rhoout_r, rhoz_r
       ! potentials from input density
-      TYPE(pw_type), POINTER                         :: vh_rspace
+      TYPE(pw_type)                                    :: vh_rspace
       TYPE(pw_p_type), DIMENSION(:), POINTER           :: vxc_rspace, vtau_rspace, vadmm_rspace
       ! efield
       TYPE(efield_berry_type), POINTER                 :: efield => NULL()
@@ -169,9 +169,8 @@ SUBROUTINE ec_env_release(ec_env)
             DEALLOCATE (ec_env%p_env)
          END IF
          ! potential
-         IF (ASSOCIATED(ec_env%vh_rspace)) THEN
+         IF (ASSOCIATED(ec_env%vh_rspace%pw_grid)) THEN
             CALL pw_release(ec_env%vh_rspace)
-            DEALLOCATE (ec_env%vh_rspace)
          END IF
          IF (ASSOCIATED(ec_env%vxc_rspace)) THEN
             DO iab = 1, SIZE(ec_env%vxc_rspace)

diff --git a/src/ec_environment.F b/src/ec_environment.F
@@ -147,7 +147,6 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section)
       NULLIFY (ec_env%matrix_wz)
       NULLIFY (ec_env%z_admm)
       NULLIFY (ec_env%p_env)
-      NULLIFY (ec_env%vh_rspace)
       NULLIFY (ec_env%vxc_rspace)
       NULLIFY (ec_env%vtau_rspace)
       NULLIFY (ec_env%vadmm_rspace)

diff --git a/src/energy_corrections.F b/src/energy_corrections.F
@@ -152,6 +152,7 @@ MODULE energy_corrections
                                               qs_kind_type
    USE qs_kinetic,                      ONLY: build_kinetic_matrix
    USE qs_ks_methods,                   ONLY: calc_rho_tot_gspace
+   USE qs_ks_reference,                 ONLY: ks_ref_potential
    USE qs_ks_types,                     ONLY: qs_ks_env_type
    USE qs_linres_kernel,                ONLY: hfx_matrix
    USE qs_mo_methods,                   ONLY: calculate_subspace_eigenvalues,&
@@ -171,8 +172,7 @@ MODULE energy_corrections
    USE qs_rho_types,                    ONLY: qs_rho_get,&
                                               qs_rho_type
    USE qs_vxc,                          ONLY: qs_vxc_create
-   USE response_solver,                 ONLY: ks_ref_potential,&
-                                              response_calculation,&
+   USE response_solver,                 ONLY: response_calculation,&
                                               response_force
    USE string_utilities,                ONLY: uppercase
    USE task_list_methods,               ONLY: generate_qs_task_list
@@ -332,7 +332,6 @@ SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr)
 
       IF (ec_env%should_update) THEN
          CALL ec_build_neighborlist(qs_env, ec_env)
-         IF (.NOT. ASSOCIATED(ec_env%vh_rspace)) ALLOCATE (ec_env%vh_rspace)
          CALL ks_ref_potential(qs_env, &
                                ec_env%vh_rspace, &
                                ec_env%vxc_rspace, &

diff --git a/src/exstates_types.F b/src/exstates_types.F
@@ -22,6 +22,8 @@ MODULE exstates_types
    USE pw_types,                        ONLY: pw_p_type,&
                                               pw_release,&
                                               pw_type
+   USE qs_local_rho_types,              ONLY: local_rho_set_release,&
+                                              local_rho_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -31,6 +33,7 @@ MODULE exstates_types
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'exstates_types'
 
    PUBLIC :: excited_energy_type, exstate_release, exstate_create
+   PUBLIC :: exstate_matrix_release, exstate_potential_release
 
 ! *****************************************************************************
 !> \brief Contains information on the excited states energy
@@ -52,11 +55,13 @@ MODULE exstates_types
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_px1_asymm => NULL()
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_px1_admm_asymm => NULL()
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_wx1 => NULL()
-      TYPE(pw_type), POINTER                           :: vh_rspace => NULL()
+      TYPE(pw_type)                                      :: vh_rspace
       TYPE(pw_p_type), DIMENSION(:), POINTER             :: vxc_rspace => NULL()
       TYPE(pw_p_type), DIMENSION(:), POINTER             :: vtau_rspace => NULL()
       TYPE(pw_p_type), DIMENSION(:), POINTER             :: vadmm_rspace => NULL()
-      REAL(KIND=dp), DIMENSION(:, :), POINTER             :: fingerprint
+      TYPE(local_rho_type), POINTER                      :: local_rho_set => NULL()
+      TYPE(local_rho_type), POINTER                      :: local_rho_set_admm => NULL()
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: fingerprint
    END TYPE excited_energy_type
 
 CONTAINS
@@ -71,7 +76,7 @@ SUBROUTINE exstate_release(ex_env)
       CHARACTER(LEN=*), PARAMETER :: routineN = 'exstate_release', &
          routineP = moduleN//':'//routineN
 
-      INTEGER                                            :: iab, is
+      INTEGER                                            :: is
 
       IF (ASSOCIATED(ex_env)) THEN
 
@@ -89,6 +94,30 @@ SUBROUTINE exstate_release(ex_env)
             END DO
             DEALLOCATE (ex_env%cpmos)
          END IF
+         !
+         CALL exstate_matrix_release(ex_env)
+         !
+         CALL exstate_potential_release(ex_env)
+         !
+         IF (ASSOCIATED(ex_env%fingerprint)) THEN
+            DEALLOCATE (ex_env%fingerprint)
+            NULLIFY (ex_env%fingerprint)
+         END IF
+
+         DEALLOCATE (ex_env)
+
+      END IF
+
+   END SUBROUTINE exstate_release
+
+! **************************************************************************************************
+!> \brief ...
+!> \param ex_env ...
+! **************************************************************************************************
+   SUBROUTINE exstate_matrix_release(ex_env)
+      TYPE(excited_energy_type), POINTER                 :: ex_env
+
+      IF (ASSOCIATED(ex_env)) THEN
          IF (ASSOCIATED(ex_env%matrix_pe)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_pe)
          NULLIFY (ex_env%matrix_pe)
          IF (ASSOCIATED(ex_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_hz)
@@ -105,42 +134,58 @@ SUBROUTINE exstate_release(ex_env)
          NULLIFY (ex_env%matrix_px1_admm_asymm)
          IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
          NULLIFY (ex_env%matrix_wx1)
-         !
-         IF (ASSOCIATED(ex_env%vh_rspace)) THEN
+      END IF
+
+   END SUBROUTINE exstate_matrix_release
+
+! **************************************************************************************************
+!> \brief ...
+!> \param ex_env ...
+! **************************************************************************************************
+   SUBROUTINE exstate_potential_release(ex_env)
+      TYPE(excited_energy_type), POINTER                 :: ex_env
+
+      INTEGER                                            :: iab
+
+      IF (ASSOCIATED(ex_env)) THEN
+         IF (ASSOCIATED(ex_env%vh_rspace%pw_grid)) THEN
             CALL pw_release(ex_env%vh_rspace)
-            DEALLOCATE (ex_env%vh_rspace)
          END IF
          IF (ASSOCIATED(ex_env%vxc_rspace)) THEN
             DO iab = 1, SIZE(ex_env%vxc_rspace)
                CALL pw_release(ex_env%vxc_rspace(iab)%pw)
                DEALLOCATE (ex_env%vxc_rspace(iab)%pw)
             END DO
             DEALLOCATE (ex_env%vxc_rspace)
+            NULLIFY (ex_env%vxc_rspace)
          END IF
          IF (ASSOCIATED(ex_env%vtau_rspace)) THEN
             DO iab = 1, SIZE(ex_env%vtau_rspace)
                CALL pw_release(ex_env%vtau_rspace(iab)%pw)
                DEALLOCATE (ex_env%vtau_rspace(iab)%pw)
             END DO
             DEALLOCATE (ex_env%vtau_rspace)
+            NULLIFY (ex_env%vtau_rspace)
          END IF
          IF (ASSOCIATED(ex_env%vadmm_rspace)) THEN
             DO iab = 1, SIZE(ex_env%vadmm_rspace)
                CALL pw_release(ex_env%vadmm_rspace(iab)%pw)
                DEALLOCATE (ex_env%vadmm_rspace(iab)%pw)
             END DO
             DEALLOCATE (ex_env%vadmm_rspace)
+            NULLIFY (ex_env%vadmm_rspace)
          END IF
-         IF (ASSOCIATED(ex_env%fingerprint)) THEN
-            DEALLOCATE (ex_env%fingerprint)
-            NULLIFY (ex_env%fingerprint)
+         IF (ASSOCIATED(ex_env%local_rho_set)) THEN
+            CALL local_rho_set_release(ex_env%local_rho_set)
+            NULLIFY (ex_env%local_rho_set)
+         END IF
+         IF (ASSOCIATED(ex_env%local_rho_set_admm)) THEN
+            CALL local_rho_set_release(ex_env%local_rho_set_admm)
+            NULLIFY (ex_env%local_rho_set_admm)
          END IF
-
-         DEALLOCATE (ex_env)
-
       END IF
 
-   END SUBROUTINE exstate_release
+   END SUBROUTINE exstate_potential_release
 
 ! **************************************************************************************************
 !> \brief Allocates and intitializes exstate_env

diff --git a/src/hfx_admm_utils.F b/src/hfx_admm_utils.F
@@ -64,8 +64,8 @@ MODULE hfx_admm_utils
         do_admm_aux_exch_func_bee, do_admm_aux_exch_func_bee_libxc, do_admm_aux_exch_func_default, &
         do_admm_aux_exch_func_default_libxc, do_admm_aux_exch_func_none, &
         do_admm_aux_exch_func_opt, do_admm_aux_exch_func_opt_libxc, do_admm_aux_exch_func_pbex, &
-        do_admm_aux_exch_func_pbex_libxc, do_potential_coulomb, do_potential_long, &
-        do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_short, &
+        do_admm_aux_exch_func_pbex_libxc, do_admm_aux_exch_func_sx_libxc, do_potential_coulomb, &
+        do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_short, &
         do_potential_truncated, xc_funct_no_shortcut, xc_none
    USE input_section_types,             ONLY: section_vals_duplicate,&
                                               section_vals_get,&
@@ -1899,6 +1899,7 @@ SUBROUTINE create_admm_xc_section(x_data, xc_section, admm_env)
 
       ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc .OR. &
                admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc .OR. &
+               admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc .OR. &
                admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
 #if defined(__LIBXC)
          IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc) THEN
@@ -1907,6 +1908,8 @@ SUBROUTINE create_admm_xc_section(x_data, xc_section, admm_env)
             name_x_func = 'GGA_X_OPTX'
          ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
             name_x_func = 'GGA_X_B88'
+         ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc) THEN
+            name_x_func = 'LDA_X'
          END IF
          !primary basis
          CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", &

diff --git a/src/input_constants.F b/src/input_constants.F
@@ -850,7 +850,8 @@ MODULE input_constants
                                                do_admm_aux_exch_func_default_libxc = 5, &
                                                do_admm_aux_exch_func_pbex_libxc = 6, &
                                                do_admm_aux_exch_func_opt_libxc = 7, &
-                                               do_admm_aux_exch_func_bee_libxc = 8
+                                               do_admm_aux_exch_func_bee_libxc = 8, &
+                                               do_admm_aux_exch_func_sx_libxc = 9
 
    ! Define scaling of the exchange correction
    INTEGER, PARAMETER, PUBLIC               :: do_admm_exch_scaling_none = 0, &